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EXECUTIVE  SUMMARY 


Microsystems  or  micro-electro-mechanical  systems  (MEMS)  are  an  emerging  technology  having 
potential  applications  in  military,  environmental,  medical  and  industrial  applications.  Though  the 
concept  of  MEMS  has  been  known  for  more  than  20  years,  the  commercialization  of  MEMS 
technology  has  not  progressed  as  fast  as  IC  chip  technology.  The  lack  of  understanding  of 
physical  phenomena  and  their  interactions  in  microsystems  has  been  a  major  technical  barrier  in 
the  use  of  MEMS  devices.  Most  microdevices  rely  on  thermal,  mechanical,  fluid,  electric,  and 
magnetic  interactions  for  performing  the  intended  function.  Recent  advances  have  been  made  in 
understanding  each  of  these  effects  individually.  However,  for  most  devices,  the  coupled  effects 
of  these  phenomena  are  not  well  understood. 

Although  the  viability  of  MEMS  devices  for  various  applications  has  been  demonstrated  in 
research  laboratories,  many  of  these  devices  have  yet  to  become  commercial  products.  The 
current  procedure  of  trial-and-error  testing  and  fabrication  of  devices  makes  them  very  expensive 
to  design.  In  many  instances  the  performance  and  reliability  are  far  from  optimal.  The 
availability  of  advanced  CAD  tools  which  could  aid  microsystem  designers  and  manufacturers  in 
analyzing/designing  microsystems  would  result  in  a  significant  reduction  in  the  extent  of 
physical  testing  that  needs  to  be  done  in  order  to  prototype  a  device.  Advanced  CAD  tools  are  a 
key  to  higher  performance/reliability,  reduced  costs,  shorter  prototyping  cycles  and  improved 
time-to-market. 

The  study  described  in  this  report  focuses  on  the  development  of  a  CAD  tool,  CFD- 
ACE+MEMS,  for  simulating  microsystem  operation.  CFD-ACE+MEMS  has  models  to  simulate 
coupled  effects  of  flow,  thermal,  mechanical,  electric,  and  magnetic  phenomena.  The  models  are 
accessible  through  geometry,  visualization,  and  problem  set-up  software.  The  simulation  abilities 
of  the  CAD  tool  have  been  demonstrated  on  many  well-characterized  microsystems  and  test 
cases.  Also,  through  experimental  measurements  conducted  at  Samoff  Corporation  (under  a  sub¬ 
contract),  the  CAD  tool  has  been  verified  using  measured  flow  field  and  mechanical 
displacements. 

CFD-ACE+MEMS  is  based  on  a  flow,  thermal,  mechanical  CAD  tool,  CFD-ACE+  which  has 
been  modified  to  include  electric  and  magnetic  physical  models.  The  structural  mechanics  model 
is  now  supplied  with  additional  forces  such  as  the  electrostatic  pressure  forces,  magnetic  moment 
forces,  and  the  Lorentz  force  on  a  conducting  body.  The  heat  model  now  has  additional  heat 
sources  due  to  joule  heating  both  from  conduction  currents  in  resistive  materials  and  inductive 
heating  due  to  time  varying  magnetic/  electric  fields.  The  flow  module  has  an  additional  Lorentz 
force  on  flows  of  (electrically)  conducting  fluids.  The  volume  of  fluids  (VOF)  model  is  supplied 
with  an  electrostatic  force  that  opposes  the  surface  tension  on  the  free  surface.  The  electric  model 
also  calculates  source  currents  for  use  by  the  magnetic  model. 


The  electric  model  solves  Poisson’s  equation  or  a  conduction  equation  for  the  electric  potential. 
From  the  electric  potential  the  electric  field,  capacitance,  electrostatic  pressure  forces,  and  the 
conduction  currents  are  calculated.  There  are  two  computational  methods  available  for  solving 
the  electric  equations.  One  is  the  finite  volume  method  (FVM)  and  the  other  is  the  boundary 
element  method  (BEM).  The  FVM  uses  finite  difference  and  the  Divergence  Theorem  to 


vii 


discretize  the  solution  space.  The  resulting  difference  equation  is  solved  using  a  linear  solver 
such  as  an  iterative  conjugate  gradient  solver  (CGS).  The  BEM  solves  the  integral  form  of  the 
electric  equations  at  boundary  surfaces.  The  boundary  surface  integrals  are  approximated  using 
techniques  similar  to  the  multipole  method.  The  solution  to  the  resulting  linear  set  of  equations  is 
obtained  using  a  solver  based  on  the  generalized  minimal  residual  (GMRES)  method. 

The  FVM  requires  a  volume  mesh  and  calculates  the  electric  potential  (and  the  electric  field)  at 
every  cell  (volume)  in  the  mesh.  The  BEM  requires  only  a  surface  mesh  (although  volume 
information  is  required  in  defining  the  different  “domains”)  and  calculates  the  electric  potential 
(and  the  electric  field)  only  at  boundary  and  interface  faces.  The  BEM  is  very  naturally  applied 
to  unbounded  problems  unlike  the  FVM  method.  The  BEM  method  has  higher  computational 
cost  (speed,  memory)  per  boundary  element  than  FVM  does  per  volume  element.  This  does  not 
mean  BEM  is  slower  since  many  problems  require  much  fewer  BEM  face  elements  than  FVM 
volume  elements,  where  the  BEM  is  faster  than  the  FVM. 

The  magnetic  model  solves  a  vector  laplacian  equation  for  the  magnetic  vector  potential  either  in 
time  or  assuming  sinusoidal  steady  state  (frequency  domain).  From  the  magnetic  vector 
potential,  the  magnetic  field,  electric  field,  Lorentz  force,  torques  on  magnetic  moments,  and 
eddy  currents  are  calculated.  The  magnetic  model  equations  are  solved  using  only  the  FVM. 

All  of  the  physical  models  (fluid,  heat,  structural,  electric,  magnetic)  are  coupled  either  (i) 
loosely,  or  (ii)  implicitly.  The  loose  coupling  is  through  boundary  and  volume  conditions  being 
passed  between  models  during  sequential  execution.  Each  physical  model  is  solved  separately 
from  the  others  using  only  boundary  conditions  and  volume  conditions  the  other  models 
calculated  during  their  previous  solution.  When  passing  the  information  between  models,  linear 
relaxation  may  be  applied  to  the  relevant  parameters  to  improve  convergence.  While  this  linking 
is  efficient  for  many  systems,  it  will  fail  when  the  coupling  becomes  too  stiff.  This  occurs  quite 
often  in  solid/fluid  systems  with  small  clearances,  i.e.,  elastohydrodynamics.  Coupled 
electric/structural  systems  with  small  clearances  will  also  exhibit  tight  coupling.  In  these 
instances,  a  certain  degree  of  implicit  coupling  between  the  models  is  necessary  to  obtain  a 
stable,  converged  solution. 

One  advanced  implicit  coupling  methodology  the  “multi-level  Newton  method”  being  developed 
by  Jacob  White’s  research  group  at  MIT  (Aluru  and  White,  1997)  has  been  implemented  to 
ensure  stable  and  robust  convergence.  The  method  is  general  enough  so  that  it  is  applicable  to 
any  physical  model  (structures,  fluids,  heat,  electromagnetics)  and  any  solution  method  (FVM, 
BEM,  FEM).  The  technique  is  a  black  box  approach  in  that  the  models  to  be  linked  need  only 
have  a  clear  set  of  input  and  output  values  upon  which  they  mutually  depend.  No  complicated 
links  need  to  be  established  between  the  models.  The  multi-level  Newton  method 
implementation  in  CFD-ACE+MEMS  has  been  demonstrated  for  an  electric/structural  coupled 

problem. 

Fundamental  validation  of  the  models  was  carried  out  by  comparing  the  CAD  tool  calculations 
with  analytical  solutions  for  simple  benchmark  cases  and  published  data  of  operating 
microsystems  in  the  literature.  Also,  a  series  of  parametric  studies  have  been  performed  to 
evaluate  the  model  characteristics  (numerical  stability  and  physical  accuracy)  over  a  range  of 
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conditions.  The  types  of  studies  carried  out  using  the  developed  CAD  tool  were  diverse.  Some 
of  these  are  listed  below: 

□  The  BEM  was  tested  for  heat  conduction  inside  a  sphere,  cylindrical  vessel  deformation 
under  surface  pressure,  and  compared  against  NASTRAN's  boundary  element  method 
package. 

□  Parallel  plates,  concentric  spheres,  and  sphere  of  space  charge  benchmarks  were  used  to 
compare  FVM  and  BEM  performance. 

□  Electrostatic  loading  of  mechanical  structures  was  modeled  for  microsystems  such  as  a 
doubly  clamped  beam,  an  accelerometer,  a  high  frequency  resonator,  an  electrostatic 
torsional  micromirror,  a  micromotor,  a  linear  lateral  resonator  comb  drive,  an  angular 
resonator  comb  drive,  and  a  fluid  damped  beam. 

□  Electrostatic  extraction  of  a  conductive  fluid  from  bath  was  calculated. 

□  The  magnetic  field  due  to  a  circular  current  carrying  wire  was  modeled.  Shielding  of 
magnetic  field  generated  from  a  bus  bar  due  to  a  high  permeability  material  was  modeled. 
The  magnetic  field  due  to  circular  and  square  planar  coils  was  calculated  using  source 
currents  from  the  electric  model.  The  magnetic  actuation  of  a  beam  with  an  attached  piece  of 
high  permeability  material  was  calculated.  The  effect  of  a  magnetic  field  on  the  damping  of  a 
buoyancy-driven  flow  of  a  conductive  fluid  was  analyzed. 

A  communication  protocol  to  use  microtomography  experimental  data  was  developed.  This 
protocol  developed  in  collaboration  with  Samoff  Corporation  facilitates  the  transfer/conversion 
of  data  files  containing  geometry  and  flow  field  information.  The  imaged  geometry  information 
can  now  be  directly  accessed  by  CFD  Research  Corporation’s  (CFDRC’s)  models  for  performing 
the  simulations.  The  protocol  also  facilitates  the  transfer  of  information  on  velocity  fields  which 
were  used  in  a  validation  study.  The  significant  outcomes  of  Samoff  s  efforts  include:  (i)  a 
progressively  refined  methodology  for  accurate  extraction  of  device  structure  from 
micro  tomography  data;  and  (ii)  a  tool  to  perform  good  quality  surface  triangulation  of  the  device. 

A  validation  study  was  carried  out  using  flow  data,  mechanical  motion,  and  geometry  drawings 
supplied  by  Samoff  Corporation.  The  measured  flow  fields  for  three  flow  channels  were  used  to 
verify  CFD-ACE  calculation  of  the  flow  fields.  Also,  previously  imaged  mechanical  motion  of 
moving  parts  in  a  comb  drive  in  conjunction  with  a  “CIF”  file  was  used  to  verify  the 
electric/mechanical  model. 

This  project  has  also  enabled  CFDRC  to  perform  several  recent  modeling  studies  (of  industrial 
microdevices)  for  MEMS  companies  using  the  coupled  simulation  capability.  These  are 
summarized  below: 

□  Micropump  for  Honeywell:  CFDRC  performed  3-D  simulations  of  coupled  fluid-structural- 
electrostatic  phenomena  in  a  micropump  currently  being  prototyped  by  Honeywell.  These 
simulations  were  successful  in  highlighting  some  of  the  complexities  of  the  interaction 
between  the  above  phenomena  and  enabled  Honeywell  designers  to  improve  the  performance 
of  the  micropump. 
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□  Electrostatically  Activated  Beam  for  Honeywell:  Simulations  were  performed  for 
Honeywell  to  analyze  the  deflection  of  a  cantilever  beam  for  different  actuation  voltages. 
Non-linear  effects  were  also  simulated.  Mechanical  forces  (as  a  function  of  gap  height)  for 
different  voltages  were  computed. 

□  Microvalve  for  Xerox:  CFDRC  performed  simulations  of  coupled  fluid-structural- 
electrostatic  interactions  in  Xerox’s  microvalve.  The  simulations  clearly  showed  the  highly 
non-linear  nature  of  the  valve  operation  and  also  gave  insights  as  to  the  extent  of  voltage 
drop  necessary  for  optimal  performance  of  the  valve. 

□  Electrophoresis  Applications:  CFDRC  is  currently  working  with  Oak  Ridge  National  Lab., 
Caliper  Technologies,  Aclara  Biosciences,  etc.  to  demonstrate  CFDRC’s  electrophoresis 
model  that  accounts  for  interactions  between  charged  species  transport  and  the  applied 
electrostatic  field. 

CFDRC  software  is  also  in  use  by  many  MEMS  groups  such  as  Redwood  Microsystems,  Lucas 
Nova  Sensor,  Motorola,  Honeywell,  Mesoscale  Systems,  YSI,  Stanford  University,  University  of 
California  (Berkeley),  University  of  Washington,  Oak  Ridge  National  Laboratory,  etc. 

The  thermo-fluidic-mechanical-electromagnetic  model  developed  and  validated  in  this  project 
has  been  incorporated  into  a  software  environment  (consisting  of  geometry  modelers,  grid 
generators,  advanced  visualization  software.  Graphical  User  Interfaces,  etc.)  and  is  being 
commercialized  as  a  design  tool.  In  addition  to  MEMS  companies,  several  microelectronics 
companies  such  as  Wacker-Siltronic,  Aixtron,  Motorola,  Applied  Materials  and  Novellus  have 
purchased  software  licenses  from  CFDRC  for  the  coupled  fluid-thermal-structural-electric- 
magnetic  software.  Their  applications  include  electromagnetic  induction  heating,  Lorentz  force 
stabilization  of  melt  flow  in  crystal  growth,  low-pressure  plasma  transport  and  electroplating. 
CFDRC  expects  a  large  market  for  this  coupled  simulation  capability.  CFDRC’s  software  is  one- 
of-a-kind  in  its  ability  to  closely  integrate  these  models. 
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INTRODUCTION 


This  is  the  final  report  documenting  the  work  performed  during  a  two  year  DARPA/AFRL 
Project  (Contract  #  F30602-97-2-0196)  entitled  “An  Integrated  CAD  Tool  for  Thermo-Fluidic- 
Mechanical-Electrostatic  Design  of  MEMS  Devices”.  The  overall  objective  is  to  develop  a  state- 
of-the-art  integrated  Computer-Aided-Design  (CAD)  tool  for  the  analysis  and  design  of 
MicroElectroMechanical  Systems  (MEMS).  Coupled  flow,  thermal,  mechanical,  electric,  and 
magnetic  models  have  been  developed.  The  coupled  models  have  been  integrated  with  geometric 
and  visualization  tools,  design  database  tools,  and  a  graphical  user  interface  (GUI)  for  ease-of- 
use  of  the  software  which  is  being  supplied  to  MEMS  vendors.  The  CAD  tool  development  has 
been  completed  and  modeling  studies  have  been  performed  on  benchmark  cases  for  testing  and 
validation  as  well  as  on  commercial  MEMS  devices.  The  advanced  3-D  imaging  and 
visualization  capabilities  of  Samoff  Corporation  were  used  to  validate  the  software.  The  CAD 
package  has  been  evaluated  on  a  range  of  microdevices  for  use  as  a  design  and  prototyping  tool. 
The  CAD  tool  is  currently  being  commercialized  to  MEMS  manufacturers  and  vendors  and 
research  groups  in  national  laboratories  and  academia. 

1.1  Microsystems  or  MEMS  (Micro-Electro-Mechanical  Systems) 

MEMS  is  an  emerging  technology  having  potential  applications  in  military,  environmental, 
medical  and  industrial  applications.  Recent  advances  in  the  development  of  micromachined 
components  such  as  sensors,  flow  channels,  valves  and  diaphragm  pumps  have  created 
opportunities  for  application  of  MEMS  technology  in  diverse  areas.  Using  a  combination  of 
these  components,  applications  in  drug  delivery  systems,  chemical  process  control, 
environmental  detection  and  industrial  feed-back  control  systems  are  becoming  possible.  Though 
the  concept  of  MEMS  has  been  known  for  more  than  20  years,  the  commercialization  of  MEMS 
technology  has  not  progressed  as  fast  as  IC  chip  technology.  The  lack  of  understanding  of 
physical  phenomena  and  their  interactions  in  microsystems  has  been  a  major  technical  barrier  in 
the  use  of  MEMS  devices.  For  example,  a  thorough  understanding  of  fluid  dynamical  behavior 
and  static  and  dynamic  interaction  of  various  components  in  complex  systems  is  lacking.  Though 
there  are  similarities  between  MEMS  components  and  their  macroscopic  counterparts,  there  are 
several  peculiarities  that  need  to  be  identified  and  understood.  For  example,  in  recent 
experiments,  it  has  been  observed  that  the  viscosity  of  liquid  in  micro-channels  is  a  function  of 
channel  size.  Also,  increased  flow  rates  and  decreased  pressure  drops  have  been  reported  which 
are  attributed  to  “slip  walls”  at  low  Knudsen  number  flows.  Most  microdevices  rely  on  thermal, 
mechanical,  electrostatic  and  fluidic  interactions  for  performing  the  intended  function.  Recent 
advances  have  been  made  in  understanding  each  of  these  effects  individually.  However,  for  most 
devices,  the  coupled  effects  of  these  phenomena  are  not  well  understood. 

Although  the  viability  of  MEMS  devices  for  various  applications  has  been  demonstrated  in 
research  laboratories,  many  of  these  devices  have  yet  to  become  commercial  products.  In  order 
for  MEMS  devices  to  be  mass-produced  commercially,  the  technology  has  to  progress  to  a  stage 
where  the  following  criteria  are  satisfied: 

1.  The  device  should  be  efficient,  i.e.,  it  should  consume  minimum  power  to  perform  the 
desired  operation. 
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2.  The  device  should  have  good  performance  characteristics  over  a  range  of  operating 
conditions. 

3.  The  reliability  of  the  device,  i.e.,  it  has  to  perform  the  desired  operation  without  the 
occurrence  of  mechanical/structural  failure  over  the  designed  lifetime  of  the  device. 

4.  The  device  design  should  be  reasonably  simple  to  facilitate  mass  production.  This  will 
lower  the  costs  of  incorporating  MEMS  technology  to  various  systems  and  thus  widen 
the  range  of  its  applicability. 


The  current  procedure  of  trial-and-error  testing  and  fabrication  of  devices  makes  them  very 
expensive.  In  many  instances  the  performance  and  reliability  are  far  from  optimal.  The  use  of 
advanced  CAD  tools  will  result  in  significant  reduction  in  the  extent  (and  hence  the  cost)  of 
physical  testing  that  needs  to  be  done  in  order  to  prototype  a  device.  Advanced  CAD  tools  will 
be  the  key  to  higher  performance/reliability,  reduced  costs,  shorter  prototyping  cycles  and 
improved  time-to-market. 

1.2  Importance  of  a  MEMS  CAD  Tool 

Although  the  viability  of  MEMS  devices  for  various  applications  has  been  demonstrated  in 
research  laboratories,  many  of  these  devices  have  yet  to  become  commercial  products.  Many 
design  issues  relating  to  device  performance/reliability,  mass  production,  packaging,  etc.,  have 
not  yet  been  addressed  satisfactorily  for  many  of  the  MEMS  devices  currently  being  explored  in 
research  laboratories.  To  achieve  these  objectives,  innovative  design  ideas  need  to  be  tested  in  a 
short  period  of  time.  Physical  testing  and  prototyping  of  these  devices  tends  to  increase  the 
design  cycle  time  and  costs  associated  with  transitioning  the  technology  to  the  commercial 
sector.  In  this  regard,  advanced  CAD  tools  will  play  a  significant  role  in  complementing  physical 
testing  and  will  aid  in  the  experimentation  and  development  of  new  designs  to  ensure 
performance  and  reliability. 

As  mentioned  before,  most  microdevices  rely  on  thermal,  fluidic,  mechanical  and  electrostatic 
interactions  for  performing  their  intended  function.  Although  advanced  models  are  available  to 
analyze  and  simulate  each  of  these  interactions,  no  tool  is  currently  available  to  model  the 
coupled  effects  of  these  phenomena.  Most  existing  models  are  able  to  assess  individual  effects  of 
the  above  phenomena.  For  example,  MEMCAD  is  able  to  perform  advanced  electrostatic 
analysis  of  MEMS  devices.  Codes  such  as  FIDAP  and  CFD-ACE  are  able  to  analyze  the  fluidic 
and  thermal  effects  in  such  systems.  Software  such  as  NASTRAN  is  able  to  do  a  thorough 
structural  analysis  of  the  system.  ANSYS  now  has  a  multi-physics  capability.  However,  a  model 
that  is  able  to  analyze  the  coupled  effects  of  these  phenomena  (flow,  thermal,  structural,  electric, 
and  magnetic)  does  not  exist. 

Many  MEMS  groups  have  used  modeling  studies  to  support  the  design  of  various  devices.  Some 
of  these  studies  are  summarized  below.  Zengerle  and  Richter  (1994)  and  Gravesen  (1993)  have 
discussed  several  aspects  of  flow  simulation  in  micro-fluidic  components.  Most  of  the 
simulations  reported  in  the  literature  utilized  very  simple  0-D  or  1-D  semi-empirical  models. 
More  recent  fundamental  Navier-Stokes  based  studies  were  performed  by  Vollmer  (1994)  on 
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2-D  fluidic  amplifiers,  by  Arkilic  (1994)  (MIT)  on  1-D  gas  flow  in  channels  and  by  Mehregany 
(1993)  on  2-D  flow  in  micro-motors.  A  review  of  the  state-of-the-art  in  computational  microfluid 
dynamics  has  been  performed  by  Przekwas  (1995). 

Early  model  development  efforts  for  MEMS  include  initiatives  such  as  MEMCAD,  at  MIT 
(Senturia  et  al.,  1992  ),  CAEMEMS  at  the  University  of  Michigan  (Zhang  et  ah,  1990),  PUSI, 
MEMS  design  automation  at  CalTech.  Recent  development  efforts  (under  the  DARPA  sponsored 
Composite  CAD  program)  include  organizations  such  as  CFDRC,  Microcosm,  Coyote  Systems, 
Tanner  Research,  Analogy,  etc.  Most  efforts  have  integrated  different  analysis  tools  (from 
different  disciplines)  into  design  environments.  Appropriate  information  on  boundary  conditions 
is  communicated  between  the  models  and  the  solution  is  obtained  by  iterating  between  the 
various  models.  However,  depending  on  the  complexity  and  non-linearity  of  coupling  between 
the  various  phenomena,  this  coupling  approach  may  not  result  in  a  converged  solution  for  the 
overall  system.  In  cases  where  the  coupling  is  highly  non-linear  and  stiff  (of  varying  time 
scales),  an  explicit  linking  of  the  models  will  not  yield  a  stable  and  converged  solution.  For 
example,  the  deflection  of  a  structure  in  a  MEMS  device  will  depend  on  the  thermal,  fluid 
pressure  and  electrostatic  loads  on  the  structure.  However,  as  the  structure  begins  to  deflect,  each 
of  the  loads  will  change  as  a  result  of  the  deflection.  In  order  to  reach  the  stable  equilibrium 
position  of  the  structure  under  the  action  of  all  the  loads,  each  of  the  models  have  to  be  coupled 
in  an  implicit  manner.  Otherwise,  convergence  will  not  be  achieved  regardless  of  the  number  of 
iterations  used  between  the  models. 

Implicit  coupling  of  models  implies  the  solution  of  one  model  depends  directly  on  the  solution  of 
the  other  models.  That  is,  the  solution  is  a  function  of  the  solution  of  the  other  methods  at 
solution  points  within  the  discretized  computational  grid.  Some  examples  of  implicit  coupling 
are  solving  the  models  simultaneously  in  one  large  matrix  or  in  some  way  having  the  predictor 
values  a  function  of  the  coupled  models. 

1.3  Focus  of  Proposed  Work 

Ongoing  modeling  activities  (such  as  development  of  MEMCAD)  have  sought  to  integrate 
specialized  software  from  several  disciplines  through  appropriate  exchange  of  information  (on 
boundary  conditions).  However,  this  approach  has  major  technical  and  commercial  drawbacks. 
MEMS  processes  are  highly  non-linear  and  the  semi-explicit  coupling  between  the  disciplines 
will  not  yield  a  converged  solution.  For  these  systems,  it  is  necessary  to  solve  all  of  the  relevant 
processes  in  a  highly  coupled  implicit  manner  in  order  to  obtain  a  stable,  converged  solution. 
Such  a  capability  is  not  currently  available  and  is  a  must  for  prototyping  MEMS  devices.  Also, 
the  use  of  a  design  environment  with  loosely  coupled  software  is  not  economically  feasible.  The 
user  will  have  to  pay  the  costs  associated  with  license  fees  for  each  of  the  software  packages 
(such  as  FIDAP,  NASTRAN,  PATRAN,  etc.)  and  still  will  not  get  adequate  technical  support  for 
the  integrated  capability. 

The  goal  of  the  proposed  work  was  to  develop  a  comprehensive  multi-disciplinary  simulation 
capability  within  a  single  environment.  Depending  on  the  problem  being  solved,  this 
environment  will  provide  options  for  either  loose  coupling  or  tight  (implicit)  coupling  between 
the  different  modules.  The  advantage  of  this  approach  is  that  it  enables  the  developer  (i.e.. 
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CFDRC)  to  provide  the  software  at  low  cost  to  the  designer.  Additionally,  the  support  of  the 
software  will  be  done  by  a  single  entity. 

The  starting  point  for  this  development  effort  was  CFD-ACE+,  a  commercial  multi-disciplinary 
software  developed  by  CFDRC.  CFD-ACE+  has  capabilities  to  simulate  coupled  effects  of  fluid 
flow,  heat/mass  transfer  and  chemistry  in  conjunction  with  stress  analysis  and  structural 
deformation.  Under  this  project,  physical  models  for  electrostatics,  magnetostatics  and 
electromagnetics  were  developed  and  coupled  into  CFD-ACE+.  The  software  has  been  tested 
and  verified  on  several  MEMS  applications  involving  coupling  between  the  modules.  Specific 
examples  include  : 

□  Electrostatic-mechanical  coupling  in  accelerometers,  comb  drives,  etc. 

□  Electrostatic-mechanical-fluidic  coupling  in  microvalves. 

□  Electrostatic-fluidic  coupling  in  electrophoresis,  electrostatic  extraction  of  droplets  from  a 
fluid  bath,  etc. 

□  Magnetostatic-fluidic  coupling  in  magnetic  stabilization  of  fluid  flow. 

□  Electrostatic-thermal-fluidic-mechanical  interactions  in  a  cantilever  beam. 

□  Magnetic-mechanical  interactions  in  a  beam 

The  completion  of  this  project  has  resulted  in  a  new  software  environment,  CFD- ACE+MEMS , 
that  is  capable  of  simulating  highly  complex  physical  interactions  in  microsystems. 
CFD- ACE+MEMS  has  set  the  stage  for  further  development  in  the  areas  of  reduced  models  and 
parametric  models  that  will  enable  complete  system  analysis.  Those  efforts  are  DARPA  BAA 
97-17:  “Generation  of  Reduced  Parametric  Models  of  Microdevices  from  High  Fidelity  Tools  for 
System  Level  Composite  CAD”  and  DARPA  BAA  97-39:  “Mixed-Dimensionality  VLSI-Type 
Configurable  Simulation  Tools  for  Virtual  Prototyping  of  Biomicrofluidic  Devices  and 
Integrated  Systems.” 

1.4  Project  Tasks  and  Accomplishments 

This  section  describes  the  tasks  that  were  proposed  under  this  project  along  with  the 
corresponding  accomplishments  for  each  task.  The  two  phases  of  work  which  have  been 
completed  in  this  project  include: 

PHASE  Is  DEVELOPMENT  OF  A  THERMO-FLUIDIC-MECHANICAL-ELECTRO¬ 
MAGNETIC  MODEL  FOR  MEMS  DEVICE  OPERATION 

Task  1.  Developed  and  Adapted  Electromagnetic  Solvers  (Electric  and  Magnetic)  for 
MEMS  Applications;  There  is  now  electric  (FVM  and  BEM  methods)  and  magnetic  solvers 
(FVM  method  only)  available  in  the  CAD  environment  (CFD-ACE+MEMS)  which  are  coupled 
to  the  flow,  heat,  and  structural  physical  models.  The  electromagnetic  solvers  have  been  tested 
and  refined  to  commercial  grade. 

Task  2.  Adanted  the  Thermo-Fluidic-Mechanical  Code  to  Include  the  Electromagnetic 
Models:  The  existing  CAD  tool  CFD-ACE+  has  been  modified  to  interface  with  the  electrostatic 
and  magnetic  models.  The  structural  mechanics  model  is  now  supplied  with  additional  forces 
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such  as  the  electrostatic  pressure  forces,  magnetic  moment  forces,  and  the  Lorentz  force  on  a 
conducting  body.  The  heat  model  now  has  additional  heat  sources  due  to  joule  heating  both  from 
conduction  currents  in  resistive  materials  and  inductive  heating  due  to  time  varying 
magnetic/electric  fields.  The  flow  module  has  an  additional  Lorentz  force  on  conductive  flows. 
The  volume  of  fluids  (VOF)  model  is  supplied  with  an  electrostatic  force  which  opposes  the 
surface  tension  on  the  free  surface.  The  electrostatic  model  can  calculate  source  currents  for  use 
by  the  magnetic  model.  The  thermo-fluidic-mechanical-electromagnetic  model  developed  and 
validated  in  this  project  has  been  incorporated  into  a  software  environment  (consisting  of 
geometry  modelers,  grid  generators,  advanced  visualization  software,  Graphical  User  Interfaces, 
etc.)  and  is  being  commercialized  as  a  design  tool. 

Task  3.  Implemented  Implicit  Coupling  Between  the  Models:  A  general  implicit  coupling 
scheme,  referred  to  as  the  multi-level  Newton  method  (Aluru  and  White,  1997),  has  been 
implemented  in  the  CAD  tool  and  demonstrated  for  mechanical-electrostatic  coupling.  The 
method  is  general  enough  so  that  it  is  applicable  to  any  physical  model  (structures,  fluids,  heat, 
electromagnetics)  and  any  solution  method  (BEM,FVM,FEM). 

Task  4.  Developed  a  Communication  Protocol  with  Microtomography  Experiments:  This 
protocol  developed  in  collaboration  with  Samoff  Corporation  facilitates  the  transfer/conversion 
of  data  files  containing  geometry  information.  The  imaged  geometry  information  can  now  be 
directly  accessed  by  CFDRC  models  for  performing  the  simulations.  The  protocol  also  facilitates 
the  transfer  of  information  on  velocity  fields  which  were  used  in  the  validation  study.  The 
significant  outcomes  of  Samoff  s  efforts  are:  (i)  a  progressively  refined  methodology  for 
accurate  extraction  of  device  structure  from  microtomography  data;  and  (ii)  a  tool  to  perform 
good  quality  surface  triangulation  of  the  device. 

Task  5.  Carried  Out  Validation  and  Parametric  Studies;  Fundamental  validation  of  the 
models  was  carried  out  comparing  the  predictions  (for  simple  systems)  with  published  data  in  the 
literature.  A  series  of  parametric  studies  has  been  performed  to  evaluate  the  model  characteristics 
(numerical  stability  and  physical  accuracy)  over  a  range  of  conditions. 

PHASE  2:  A  PPL  I  CATION  AERIFICATION  OF  PHYSICAL  MODELS  ON 

INDUSTRIAL  DEVICES 

Task  1.  Selected  Microsystems  to  be  Modeled:  Selected  microsystems  to  be  modeled  based  on 
the  following  criteria:  (i)  the  device  operation  should  rely  on  physical  phenomena  that  can  be 
simulated  by  the  model,  (ii)  the  device  should  be  sufficiently  experimentally  characterized  to 
enable  modeling. 

Task  2.  Performed  Verification  Simulations  of  Selected  Microsystems;  The  microfluidic 
simulations  of  Samoff  devices  have  been  completed.  Electrostatic  model  verification  of  an 
angular  resonator  comb  drive  has  been  completed.  The  parameters  that  were  tested  are  (i) 
geometry  of  the  device,  (ii)  the  materials  used  in  the  device,  (iii)  flow  conditions,  and  (iv) 
magnitude  of  structural  deflection  of  microsystems. 
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Task  3.  Visualization  of  Device  Performance:  The  operation  of  microsystems  were  imaged 
and  visualized  at  Samoff  Corporation.  Specifically,  the  flow  channel  geometries  have  been 
imaged  and  the  velocity  fields  of  flows  in  these  channels  have  been  measured.  Previously 
imaged  mechanical  motion  of  moving  parts  in  a  comb  drive  was  used  to  verify  the 
electric/mechanical  model. 

1.5  Outline  of  the  Report 

Chapter  2  describes  the  fundamental  formulation  and  governing  equation  for  each  of  the  physical 
models  developed  and  adapted  during  the  course  of  this  project.  Chapter  3  describes  the 
application  of  the  finite  volume  method  (FVM)  and  the  boundary  element  method  (BEM)  for 
computing  electrostatic  fields.  The  integration  of  the  physical  models  into  a  commercial  CAD 
environment,  CFD- ACE+MEMS ,  is  described  in  Chapter  4.  Detailed  validation  and  model 
verification  studies  are  presented  in  Chapter  5.  Chapter  6  describes  the  application  of  CFD- 
ACE+MEMS  for  simulating  coupled  multi-disciplinary  phenomena  in  industrial  microsystems. 
The  implicit  coupling  procedure  between  the  individual  modules  is  presented  in  Chapter  7. 
Chapter  8  describes  the  imaging  and  visualization  activity  performed  at  Samoff  Corporation. 
Chapter  9  concludes  with  a  summary  of  work  accomplished,  commercialization  efforts  and  plans 
for  future  development. 
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2.  DESCRIPTION  OF  PHYSICAL  MODELS 

This  chapter  describes  the  physical  models  in  the  CAD  tool.  First,  the  electromagnetic  models 
(electric  and  magnetic  models)  developed  under  this  project  will  be  described.  During  the 
discussion,  emphasis  will  be  placed  on  how  the  electromagnetic  models  have  been  coupled  to 
existing  flow,  heat,  and  structural  models.  The  structural  mechanics  model  is  now  supplied 
additional  forces  such  as  the  electrostatic  pressure  forces,  magnetic  moment  forces,  and  the 
Lorentz  force  on  a  conducting  body.  The  heat  model  now  has  additional  heat  sources  due  to  joule 
heating  both  from  conduction  currents  in  resistive  materials  and  inductive  heating  due  to  time 
varying  magnetic/electric  fields.  The  flow  module  has  an  additional  Lorentz  force  in  the 
momentum  equation  for  conductive  flows.  The  electric  model  can  calculate  source  currents  for 
use  by  the  magnetic  model. 

2.1  Electric  Model 

The  electric  model  involves  the  solution  of  the  electric  potential.  The  model  solves  two  basic 
equations,  either  Poisson’s  equation  (Eq.  (2.1))  or  a  conduction  current  (Eq.  (2.2))  equation. 


V  •arV<j>=--^-  (2.1) 

ao 

V  •6V<j)=0  (2.2) 

where  (j>  is  the  electric  potential,  er  is  the  electric  relative  permittivity,  e0  is  the  permittivity  of 
free  space,  p  is  the  space  charge,  and  a  is  the  electric  conductivity.  From  the  electric  potential 
the  electric  field  E,  the  pressure  force,  and  the  capacitance  of  the  geometry  are  calculated.  If  the 
conduction  problem  is  solved  the  current  density  J,  virtual  force,  and  joule  heating  source  are 
also  calculated. 

Poisson’s  equation  comes  from  Gauss’  law  for  the  electric  field  (V*D  =  p),  the  constitutive 
relation  D  =  £0srE,  and  the  definition  of  the  static  electric  field  in  terms  of  the  electric  potential 
(i.e.,  E  =  -V(J)  since  VxE  =  0).  The  electric  field  which  satisfies  Poisson’s  equation  is  used  to 
calculate  pressure  forces  and  capacitances  when  an  electrostatic  assumption  is  valid  (no 
movement  of  charge). 

The  conduction  problem  is  primarily  used  to  calculate  source  currents  for  the  magnetic  model, 
joule  heating  sources  for  the  heat  module,  and  pressure  forces  for  the  volume  of  fluids  model. 
The  conduction  or  current  continuity  equation  (Eq.  (2.2))  is  a  reduced  form  of  the  continuity 
equation  (Eq.  (2.3)).  The  conduction  equation  used  here  assumes  movement  of  charge  is 
occurring  at  a  much  faster  rate  then  the  time  steps  taken  (no  time  varying  magnetic  fields, 
electric  fields,  or  space  charge),  and  that  the  movement  of  charge  is  due  only  to  conduction 
current  (no  convection  current).  Also,  implicit  in  the  conduction  equation  here  is  that  the 
electrons  have  had  enough  time  to  move  through  the  material(s)  to  a  relaxed  state  (i.e.,  steady 
state). 
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(2.3) 


V 


•J  =  - 


dn 

dt 


2.1.1  Capacitance 

Capacitance  is  a  ratio  of  separated  charge  to  electric  potential.  In  the  electric  model  a  total  charge 
is  calculated  for  every  set  of  boundaries  with  a  unique  fixed  electric  potential  <t>0.  To 

determine  the  Q  at  a  fixed  potential  boundary,  one  starts  with  the  electric  field  boundary 
conditions 


Ps 


n 


2 


1 


Figure  2-1.  Defining  Region  Numbers  and  Direction  of  Normal  for  Eq.  (2.4) 


n-(a1E1-a2E2)=ns  (2.4a) 

n  x  (Ei  -E2)=0  (2.4b) 


Eq.  (2.4a)  states  that  the  discontinuity  at  an  interface  in  the  normal  component  of  the 
displacement  flux  is  equal  to  the  surface  charge  at  the  interface.  Assume  one  region  (2)  is 
metallic  so  that  the  electric  field  goes  to  zero  there.  The  total  charge  on  a  contact  is  then  given  by 
Eq.  (2.5).  The  total  change  is  the  sum  at  each  face  with  a  fixed  potential  (<j)0)  of  the  incident 
displacement  flux  normal  to  the  fixed  electric  potential  boundary  (stjiiEni)  times  the  area  of  the 
face  (Aj). 


faces  with  be  <))0 

Qo0  = 

i  =  l 


(2.5) 


A  capacitance  is  calculated  for  each  set  of  faces  with  a  unique  potential  boundary  condition.  First 
the  total  charge  on  each  collection  of  faces  with  a  unique  fixed  potential  boundary  condition 


given  by  Eq.  (2.5)  is  calculated.  Then  for  those  face  sets  with  <j>o*0.0  the  ratio  C  = 


^ o 


<f>c 


is 


calculated. 


2.1.2  Electrostatic  Pressure  Force  (Electric/Structural) 

The  coupling  between  the  electrostatic  model  and  the  structural  mechanics  model  is  through  the 
pressure  forces.  The  electric  model  solves  Poisson’s  equation  and  uses  the  electric  fields  to 
calculate  the  pressure  force.  The  stress  model  uses  the  pressure  forces  as  boundary  conditions. 

The  force  F  exerted  by  the  total  electric  field  Etot  on  a  surface  S  with  a  surface  charge  ps  is  given 
by 
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(2.6) 


F  =  jnsEtotda  or  dF  =  nsEtotda 

s 

The  electric  field  at  the  surface  is  discontinuous  across  the  surface  interface  as  given  by  the 
electric  field  boundary  conditions  stated  in  Eq.  (2.4). 

A  good  approximation  for  the  total  electric  field  in  Eq.  (2.6)  is  to  take  the  average  of  the  electric 
field  on  either  side  of  the  interface. 


Etot*^(El+E2)  <2-7> 

Substituting  Eq.  (2.4a)  and  Eq.  (2.7)  into  Eq.  (2.6)  yields 

dF  =  nsEtotda  =  ns  —  (Ej  +  E2)da  =  n  -(a^  -  a2E2)  —  (Ej  +  E2)da  (2.8) 

So,  the  pressure  force  dF/da  exerted  by  the  total  electric  field  Etot  (approximated  as  an  electric 

field  average)  on  an  area  element  da  with  a  surface  charge  ps  separating  two  regions  (1  and  2)  is 

P  =  ^  =  n-(a1E, -a2E2)i(E, +E2)  (2.9) 

da  2 

where  the  normal  is  directed  toward  region  1  and  e  is  the  electric  permittivity. 

Eq.  (2.9)  can  be  simplified  if  one  of  the  two  regions  is  a  metal.  For  instance,  if  region  2  is  a 
metal,  no  electric  field  exists  in  the  region  (E2  =  0),  and  since  the  tangential  component  of  the 
electric  field  is  continuous  across  the  interface  (Eq.  (2.4b)),  the  tangential  component  of  the 
electric  field  in  region  1  is  also  zero  (E]t  =  0).  So,  if  region  2  is  a  metal  then  Eq.  (2.9)  reduces  to 

P  =  —  =  —  e1Einn  if  metal  in  region  2  (2.10) 

da  2 

where  Em  is  the  normal  component  of  the  electric  field  in  region  1,  £1  is  the  permittivity  in 
region  1,  and  the  normal  is  directed  toward  region  1  as  shown  in  Figure  2.1. 

2.1.3  Virtual  Force  (Electric/Volume  of  Fluids) 

A  virtual  force  calculation  was  added  to  ease  the  coupling  with  a  volume  of  fluids  (VOF)  model 
which  is  a  technique  for  resolving  free  surfaces.  Since  boundary  conditions  at  the  free  surface  are 
no  longer  directly  associated  with  a  mesh  face  a  force  calculation  that  did  not  rely  on  an 
available  list  of  boundaries  was  needed.  The  solution  was  to  calculate  the  virtual  force  at  a 
volume  due  to  the  varying  electric  energy  in  the  region.  The  electrostatic  potential  is  used  to 
calculate  a  force  at  the  surface  of  the  liquid  due  to  electrostatic  pressure  by  solving 
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Fyirt  -  ^ 


\  «N 


(2.11) 


where  s  is  the  electric  permittivity.  The  virtual  force  per  unit  volume  FVjrt  is  obtained  by  taking 
the  gradient  of  the  energy  per  unit  volume  stored  in  the  electric  field.  To  couple  electrostatics 
with  VOF  method,  VOF  is  given  access  to  electrostatic  material  properties  and  recalculates  them 
as  a  function  of  fluid  fraction.  The  conduction  problem  is  then  solved  and  the  electrostatic  virtual 
forces  at  every  cell  center  is  calculated  by  taking  the  gradient  of  the  energy  stored  in  the  electric 
field.  For  coupling  to  the  VOF  method,  the  virtual  force  is  applied  only  to  cells  which  have  a 
fractional  content  of  the  liquid.  Elsewhere  in  the  domain  the  force  is  assumed  zero.  These  virtual 
forces  are  passed  as  body  force  to  fluid  model. 

2.1.4  Electric  Current  Density  (Electric/Magnetic,  Electric/Heat) 

The  solution  of  the  conduction  equation  (Eq.  (2.2))  supplies  links  to  magnetics,  heat,  and  VOF 
models.  The  electric  potential  solution  from  the  conduction  problem  is  used  to  calculate  source 
currents  for  the  magnetics  model,  joule  heating  sources  for  the  heat  module,  and  body  forces  for 
the  volume  of  fluids  (VOF)  model. 

The  conduction  equation  provides  source  currents  (Eq.  (2.12))  for  the  magnetics  model.  This  is  a 
useful  method  for  describing  source  currents  in  that  the  user  does  not  have  to  specify  a  direction 
of  the  current,  it  enables  the  use  of  source  currents  of  arbitrary  geometry,  and  it  allows  the 
current  density  direction  to  reflect  any  changes  due  to  a  moving  grid. 


J=aE  (2-12) 

When  the  conduction  problem  is  solved  a  source  term  is  also  added  to  the  heat  equation  (J*E) 
due  to  current  flow  through  conductive  materials. 

2.1.5  Boundary  Element  Method  (BEM)  and  Finite  Volume  Method  (FVM) 

Two  different  techniques  are  available  for  solving  the  electric  equations  (Poisson  or  conduction 
equations).  One  is  the  finite  volume  method  (FVM)  and  the  other  is  the  boundary  element 
method  (BEM).  The  FVM  requires  a  volume  mesh  and  calculates  the  electric  potential  (and  the 
electric  field)  at  every  cell  (volume)  in  the  mesh.  The  BEM  requires  only  a  surface  mesh 
(although  volume  information  is  required  in  defining  the  different  “domains  )  and  calculates  the 
electric  potential  (and  the  electric  field)  only  at  boundary  and  interface  faces.  A  detailed 
discussion  of  BEM  and  FVM  is  given  in  the  next  chapter 

2.2  Magnetic  Model 

A  magnetic  model  was  developed  to  compute  forces  on  conductive  flows  and  model 
microsystems  which  rely  in  part  on  magnetic  actuation  for  operation.  The  magnetic  model  solves 
a  reduced  set  of  Maxwell’s  equations  for  a  magnetic  vector  potential.  The  magnetic  vector 
potential  is  used  to  calculate  the  magnetic  field,  forces  on  high  permeability  materials,  and 
Lorentz  forces  on  current  carrying  wires  and  conductive  flows.  The  magnetics  model  has  been 
implemented  in  2-D,  2-D  axisymmetric,  and  3-D  using  the  FVM. 
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The  magnetic  model  solves  Maxwell’s  equations  for  the  magnetic  vector  potential  A 

V2  —  =  -  Js  -  aE  +  ct  (2.13) 

p  ot 

assuming  negligible  displacement  current  —■  =  0,  isotropic  permeability  p,  and  the  Coulomb 
gauge  V •  A— 0.  The  current  sources  can  be  user  specified  (Js),  calculated  from  the  solution  of  a 

A 

conduction  problem  (ctE),  or  eddy  ( — — )  currents.  Eq.  (2.13)  is  solved  on  a  2-D  (cartesian, 

3t 

cylindrical)  or  3-D  unstructured  mesh  using  the  FVM.  In  2-D  the  source  currents  are  scalars  with 
assumed  direction,  Js  =  Jzz  (cartesian),  or  Js  =  J^<j)  (cylindrical).  From  the  magnetic  vector 

potential  the  magnetic  field  B  can  be  calculated  (B  =  V xA). 

For  sinusoidal  current  sources  (ac  case)  the  equation  can  be  solved  in  phasor  form  to  obtain  the 
sinusoidal  steady  state  solution 


V2  —  =  -  Js  +  jatoA  (2.14) 

The  coupling  between  the  magnetic  model  and  the  structural  mechanics  model  is  through  volume 
forces.  The  magnetic  model  solves  Eq.  (2.13)  or  Eq.  (2.14)  for  the  magnetic  vector  potential  A 
which  is  used  to  calculate  the  magnetic  field  (B  =  VxA).  The  magnetic  field  is  used  to  calculate  a 
volume  force.  The  stress  model  uses  the  volume  forces  as  body  force  conditions.  Different 
regions  such  as  current  carrying  solids/fluids  and  high  permeability  solids  require  different  force 
calculations. 

2.2.1  Magnetic  Force  on  Current  Carrying  Wire  (Magnetic/Structure) 

For  regions  where  current  sources  exist  the  incremental  force  dF  due  to  a  magnetic  field  B  on  an 
incremental  volume  element  dv  with  a  current  density  J  is  calculated  using 

dF  =  Jdv  x  B  (2.15) 

2.2.2  Force  on  a  High  Permeability  Material  (Magnetic/Structure) 

For  regions  were  magnetized  material  exist  either  the  virtual  force  per  unit  volume  or  the  torque 
on  a  magnetic  moment  can  be  calculated. 

Virtual  Force 

For  regions  were  magnetized  material  exist  the  virtual  force  per  unit  volume  on  the  volume 
element  is  calculated  by  taking  the  gradient  of  the  magnetic  energy  density  of  the  volume 
element,  Fvirt  =  VUmag.  If  electric  fields  also  exist  in  the  region  the  energy  density  of  the  electric 
field  is  included  in  the  gradient  calculation.  The  instantaneous  energy  density  of  a  magnetic  field 
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is  Uma2=  -B«H.  When  the  fields  are  in  Phasor  form  Umag  =^-Re[B*H°]  where  the  energy 
represents  the  average  over  one  cycle,  and  the  symbol  represents  the  complex  conjugate. 

Force  Due  to  a  Magnetic  Moment 

A  magnetic  moment  force  calculation  has  been  added  in  addition  to  the  virtual  force  calculation. 
A  material  with  a  magnetic  moment  M  placed  in  a  magnetic  field  H  will  experience  a  force  F 
given  by 

F  =  p0(M*V)H  (2-16) 


The  magnetic  moment  M  is  currently  assumed  proportional  to  the  applied  magnetic  field  (M  = 
XmH).  So  using  the  definition  of  a  magnetic  moment  in  terms  of  B  and  H,  M  =  B/p0  -  H  and  the 
definition  of  relative  permeability,  pr  =  (1  +  Xm)  Eq.  (2.16)  becomes 


fi-0 

F  = 

B 

•  V 

V  Pry 

B 


M'oPr 


(2.17) 


which  gives  the  force  on  a  magnetizable  material  in  terms  of  the  applied  magnetic  field. 


2.2.3  Force  on  a  Conductive  Flow  (Magnetic/Flow) 

The  Lorentz  force,  J  x  B  is  included  in  the  momentum  equation  for  the  flow  of  an  electrically 
conductive  fluid  in  a  magnetic  field.  The  current  density  J  is  given  by 

J  =  a(E+u  x  B)  (2-18) 


The  Lorentz  force  can  be  expanded  as 

F  =  a(-V(j)  +  u  x  B)  x  B  =  -  aV<j>  x  B  +  a(u  •  B)B  -  a(B  •  B)u  (2.19) 

The  last  two  terms  damp  the  flow  when  the  magnetic  field  B  and  velocity  u  are  not  aligned. 

The  current  in  the  fluid  is  calculated  from  the  conservation  of  charge.  The  divergence  of  the 
current  density  is  zero,  giving 

V  J  =  V-c(-V<j>  +  u  x  B)  =  0  (2.20) 

This  equation  is  solved  for  the  electric  potential,  <j>.  The  electrostatic  model  was  modified  to 
include  the  term  V-ct(u  x  B).  The  boundary  conditions  for  the  electrostatic  module  have  been 
modified  to  allow  for  insulated  boundaries  where  the  normal  component  of  the  current  density 

vector  is  zero. 

2.2.4  Joule  Heating  bv  Induction  (Magnetic/Heat) 

A  time  varying  magnetic  field  generates  a  time  varying  electric  field  which  generates  eddy 
current  in  conductive  materials.  The  flow  of  current  through  the  conductor  generates  heat  and  is 
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called  inductive  heating.  The  inductive  heating  can  be  thought  of  as  joule  heating  (J*E)  with  the 
conductive  currents  (J  =  ctE)  generated  by  the  time  varying  field  (E  =  -V<j)  -  3A/3t).  So  assuming 
the  electric  field  is  due  only  to  the  time  vary  magnetic  field  (E  =  jcoA)  and  the  current  is  only  a 
conduction  current  (J  =  a  A)  the  time  averaged  expression  for  joule  heating  becomes 


1  2 
—  arco 

2  r 


where  ar  is  the  real  component  of  conductivity  and  co  is  the  radian  frequency  (2nf). 


(2.21) 


2.3  Flow  Model 

The  fluid  flow  model  solves  the  continuity  equation  and  the  pressure  based  Navier-Stokes 
equations  (Patankar,  1980)  assuming  that  the  fluid  is  a  continuum. 


f(pUi)+ 


dp  d  \  . 

—  + - (pu ;  =  n 

St  ^ -  J' 


m 


^4uiuj)=~ 


dx 


OP 

dX; 


+ 


+  Pf, 


(2.22) 

(2.23) 


where  p  is  the  fluid  density,  Uj  is  the  jth  Cartesian  component  of  the  instantaneous  velocity,  m  is 
the  rate  of  mass  generation  in  the  system  (typical  examples  of  mass  generation  include  surface 
reaction,  spray  vaporization  etc.),  p  is  the  static  pressure,  Xq  is  the  viscous  stress  tensor,  and  fj  is 
the  body  force  (e.g.,  gravity,  surface  tension,  etc.). 

For  very  small  flow  dimensions  and/or  very  low  pressure  gas  flows,  the  mean-free  path  of  the 
molecules  X  becomes  comparable  to  the  characteristic  flow  dimensions  L  and  the  no-slip  wall 
conditions  associated  with  continuum  assumption  are  no  longer  valid.  For  instance,  low  pressure 
gas  flows  in  micro-fluidics.  Use  of  no-slip  walls  under  these  conditions  usually  generates  wall 
shear  that  is  overestimated.  If  the  operating  conditions  are  such  that  the  Knudsen  number  (Kn  = 
XfL)  falls  in  the  “transitional”  regime  (0.01  <  Kn  <  3.0),  a  slip-wall  boundary  condition  is  used. 
The  slip-wall  formulation  includes  effects  of  local  temperature  and  diffusive  creep  in  velocity 
slip.  A  formulation  neglecting  thermal  effects  can  also  be  used  (Rohsenow  and  Choi,  1961). 
Under  low  absolute  pressures  and  gaseous  mixtures  the  gas  transport  properties  are  not  known 
and  the  Lennard  Jones  potential  model  can  be  used  to  calculate  the  transport  properties  of 
rarefied  gas/gas  mixtures  at  low  absolute  pressure.  A  detailed  description  is  given  in  Bird  et  al., 
1960. 


The  salient  features  of  the  fluid  model  are: 

•  pressure-based,  implicit  Navier-Stokes  code 

•  incompressible  and  compressible  flows 

•  turbulence  models  including  standard  k-e 

•  steady-state  or  transient  problems 

•  multi-component  diffusion 


•  multi-species  transport 

•  multi-step  gas  phase  chemical  reactions 

•  surface  reactions 

•  slip-wall  boundary  conditions 
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2.4  Heat  Model 


The  thermal  model  solves  the  total  enthalpy  form  of  the  energy  equation 


|‘(PH)  +  — — (pUjH): 
dt  ox  j 


8  ,  9p 

— (-qj)+—  + 


dXj 


dt  <3x; 


(TijUj) 


5Xj 


(Jijhi)+Sa+pfiui 


(2.24) 


where  Jy  is  the  total  (concentration-driven  +  temperature-driven)  diffusive  mass  flux  for  species 
i,  hj  represents  the  enthalpy  for  species  i,  Sa  represents  additional  sources  due  to  surface  reaction, 
radiation,  joule  heating,  and  liquid  spray,  and  qy  is  the  j-component  of  the  heat  flux.  Fourier’s 
Law  is  employed  to  model  the  heat  flux 


(2.25) 


where  K  is  the  thermal  conductivity  and  T  is  the  temperature.  The  total  enthalpy  H  is  defined  as 


H  =  h  + 


u)u) 


(2.26) 


where  h  is  the  static  enthalpy. 

The  energy  equation  includes  unsteady,  convective,  conductive,  species  energy,  viscous 
dissipation,  work,  joule  heating  and  radiation  terms.  The  conductive  term  uses  conjugate  heat 
transfer.  The  radiation  source  term  in  the  energy  equation  is  a  solution  of  a  radiative  heat  transfer 
equation  for  emitting,  absorbing,  scattering,  gray  medium  using  either  a  discrete  ordinate  method 
(Fiveland,  1984)  or  a  surface-surface  exchange  method. 

2.5  Mechanical  Model 

The  mechanics  equations  describe  the  deformation  of  a  solid  body,  and  the  resulting  strains  and 
stresses,  under  applied  loads.  These  applied  loads  may  be  surface  forces  (e.g.  pressure)  and/or 
body  forces  (e.g.  gravity). 

2.5.1  Strain-Displacement  Relations 

In  fluid  mechanics,  we  are  typically  interested  in  measuring  the  velocity  of  fluid  particles  as  they 
pass  through  a  fixed  point  in  space,  and  thus  an  Eulerian  formulation  is  used.  In  solid 
mechanics,  on  the  other  hand,  we  use  a  Lagrangian  formulation  where  a  continuum  undergoes  a 
deformation  with  the  motion  expressed  as: 


R  =  R(r,t)  (2.27) 

which  gives  the  current  location  R  of  the  particle  that  occupied  the  position  r  at  time  t  =  0.  This 
may  be  thought  of  as  a  mapping  of  the  initial  configuration  onto  the  current  configuration.  Thus, 
we  may  write 
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r  =  x,i,  (2.28) 

R  =  r  +  u  (2.29) 

Where  u  is  the  displacement  vector  of  the  particle. 

Consider  an  infinitesimal  line  segment  on  the  undeformed  and  deformed  bodies.  The  square  of 
the  lengths  of  this  line  segment  are  given  by: 

ds2  =  dr  ■  dr  =  dxjdxj  (2.30) 

ds2  =  dR-dR  =  dxjdx ;  (2.31) 

dx  :  dx ;  J 


The  difference  between  ds2  and  ds2  is  a  measure  of  strain  in  the  body,  defined  as 

ds2  -  ds2  =  2yydxidXj  (2.32) 

where  yy  is  called  the  Green-Saint-Venant  strain  tensor.  Using  Eqs.  (2.28)  through  (2.31),  the 
strain  tensor  is  written  as 


Yij 


~  (ui,j  +  uj,i  +  um,ium,j) 


(2.33) 


This  representation  is  referred  to  as  the  Lagrangian  strain  tensor  since  the  initial  coordinates  were 
used  as  the  basis. 


If  the  displacement  gradients  Ujj  are  small  in  comparison  to  unity,  the  product  term  in  Eq.  (2.33) 
is  negligible  and  the  strain  is  represented  by  its  linear  portion. 

Tii—kj  +  uy)  <2-34> 

Such  problems  are  referred  to  as  geometrically  linear. 

2.5.2  Stress-Strain  Relations 

A  body  which  has  undergone  deformation  as  a  result  of  applied  loads  will  transmit  forces 
throughout  its  interior  to  balance  the  applied  loads.  As  a  result  of  forces  being  transmitted  from 
one  portion  of  the  continuum  to  another,  an  element  with  area  AAn  on  or  within  the  body  and 
oriented  with  normal  n,  will  be  subjected  to  a  resultant  force  AFn  and  moment  AMn.  The  vector 
force  and  moment  will  not,  in  general,  be  parallel  to  n. 
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The  Cauchy  Stress  Principle  states  that  the  ratio  AFn/  AAn  tends  to  be  a  definite  limit  as  AAn 
approaches  zero,  and  also  that  the  moment  vanishes  in  the  limiting  process.  The  stress  vector  or 
traction  is  defined  as 


=  lim 
AAn  — >0 


AFn  _  dFn 
AAn  dAn 


(2.35) 


The  traction  Tj,  acting  on  a  face  normal  to  the  Cartesian  unit  vector 
terms  of  the  Cartesian  stress  tensor  a,j  as 


ij ,  may  be  expressed  in 
(2.36) 


The  coefficients  c jy  represent  components  of  the  stress  tensor.  The  subscript  and  sign 
conventions  are  as  follows: 

•  The  first  subscript  i  refers  to  the  normal  ij  denoting  the  face  on  which  Tj  acts. 

•  The  second  subscript  j  refers  to  the  direction  ij  in  which  the  stress  acts. 

•  The  normal  components  (an,  (J22,  033)  are  positive  in  tension  and  negative  in 
compression. 

•  The  shearing  components  aij  (i  ^  j)  are  positive  if  directed  in  the  positive  Xj  direction 
while  acting  on  the  face  with  the  unit  normal  +  ij ,  or  if  directed  in  the  negative  xj 

direction  while  acting  on  the  face  with  unit  normal  -  ij . 

If  the  Cartesian  stress  tensor  ay  is  known,  this  completely  defines  the  state  of  stress  at  a  point 
(Gould,  1983). 

A  relationship  between  the  stress  tensor  ay  and  strain  tensor  yy  may  be  derived  as  follows. 

A  generalized  Hooke’s  Law  formulation  may  be  derived  by  assuming  the  existance  of  an  elastic 
potential  W  which  is  invariant  with  coordinate  transformations  of  the  strain  tensor  Yy: 

W  =  W(Yij) 

Expanding  this  in  a  Taylor  Series  about  y y  =  0  gives 

W  =  W0+EjjYjj+-EjjkmyijYkm  +  —  Ejjj^npYijYkmYnp  +■•• 


(2.37) 


(2.38) 


where  Wo  is  a  constant  and 


F  -5W  F 

ij-^,  .  ’  ‘j1™ 


dy 


ij 


a2w 


jEjj 


53W 


ijkmnp 


ij^Ykm^f  np 


,etc. 


(2.39) 
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(2.40) 


Consider  SW,  a  small  change  in  energy  due  to  Syy, 

8W  =  Ojj8yij 


The  81  coefficients  Eijki  are  called  elastic  constants.  With  application  of  stress  and  strain  tensor 
symmetry,  the  number  of  independent  equations  reduces  from  nine  to  six,  and  the  number  of 
elastic  constants  reduces  to  21  (using  symmetry  of  the  Eyk|  tensor).  Further  reductions  are 
possible  for  materials  with  one  or  more  axes  of  symmetry. 

Note  that  small  strain  does  not  necessarily  imply  small  displacement  gradients  (and  thus,  a  linear 
strain-displacement  relationship).  It  is  possible  to  have  a  large  deformation,  small  strain 
problem. 

2.5.3  Governing  Equation  of  Mechanics 

The  governing  equation  of  the  mechanics  module  is  written  in  terms  of  stress,  which  may  then  be 
related  to  displacement  via  the  stress-strain  and  strain-displacement  equations.  The  equation  is 
derived  using  a  balance  of  linear  momentum  as  follows. 


17 


Consider  a  volume  V  with  surface  area  A.  The  principle  of  linear  momentum  states 

J  fdV  +  }  TdA  =  f  piidV  (2.45) 

V  A  v 

where  f  is  the  body  force,  T  is  the  surface  traction,  p  is  the  mass  density,  and  u  is  the 
displacement  vector. 

We  write  this  in  tensor  form  using 

f=fiii  (2.46a) 

T =Tiii  =crjin  jij  (2.46b) 

u=Ujij  (2.46c) 

The  Green-Gauss  theorem  is  used  to  convert  the  area  integral  of  Eq.  (2.45)  into  a  volume 
integral: 

|aijnjiidA=Jaji(jiidV  (2.47) 

A  V 

Equations  (2.45)  and  (2.47)  give 

j(cfji,j+fi-Piii)idV  =  0  (2'48) 

V 

For  this  integral  to  vanish  for  all  arbitrary  volumes  the  integrand  must  vanish,  which  yields 

ajij+fi -phi  =°  (2-49) 


This  is  the  equation  of  motion  describing  the  distribution  of  stress  in  the  body. 


The  balance  of  angular  momentum  shows  that  the  stress  tensor  is  symmetric  (Chung,  1996),  i.e. 


°ij  -  a  ji 


(2.50) 


Thus,  Eq.  (2.49)  is  often  written  as 


°ij,j  +fi  -P^i  =0 


(2.51) 


2.5.4  Multidisciplinary  Coupling 

The  mechanics  module  is  coupled  to  other  modules  of  CFD-ACE+  through  boundary  conditions 
and  source  terms.  These  are  described  briefly  below. 
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Coupling  with  Thermal  Module  Thermoelasticity 

Temperature  changes  can  be  an  important  contributor  to  stress  in  elastic  bodies.  We  assume  a 
linear  relationship  between  the  temperature  difference  from  a  datum  value  and  the  corresponding 
strain: 


T(ii)(T)=a(i)[T-To] 


(2.52) 


where  the  parentheses  in  the  subscript  implies  no  summation.  Here  Y(ii)  (T)  is  the  linear  strain  in 
coordinate  direction  i  corresponding  to  the  temperature  difference  T  -  T0,  and  a(i)  is  the 
coefficient  of  thermal  expansion  corresponding  to  that  coordinate  direction. 

With  Eq.  (2.52),  we  generalize  the  stress-strain  relations  to 

aij  =Eijkl(Ykl  +5kla(k)[T_To])  (2-53) 


Coupling  with  Electrostatics  and  Fluids 

The  fluid  and  electrostatic  modules  are  each  coupled  to  the  mechanics  module  via  surface 
pressure  terms.  These  enter  the  system  as  boundary  conditions  to  Eq.  (2.51),  where  the  stress  at 
the  boundary  is  given  by  the  expression 


CTijnj  =  Pni  (2'54) 

at  a  boundary  with  unit  normal  n.  Here,  p  is  the  pressure  from  the  fluids  or  electrostatics 
module. 

Coupling  with  Magnetics 

The  magnetic  model  is  coupled  to  the  mechanics  model  via  a  body  force.  The  magnetic  body 
force  can  occur  in  current  carrying  wires  where  the  moving  charge  experiences  a  force  as  it 
travels  through  a  magnetic  field  (Eq.  (2.15)).  A  magnetic  force  can  also  occur  on  magnetized 
material  (Eq.  (2.17))  as  the  magnetic  field  attempts  to  align  the  magnetic  moment. 

The  body  force  vector  f  is  then  given  directly  to  Eq.  (2.51)  as  a  source  term. 
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3.  BOUNDARY  ELEMENT  AND  FINITE  VOLUME  METHODS 

Two  different  techniques  are  available  for  solving  the  electric  equations  (Poisson  or  conduction 
equations).  One  is  the  finite  volume  method  (FVM)  and  the  other  is  the  boundary  element 
method  (BEM).  The  FVM  requires  a  volume  mesh  and  calculates  the  electric  potential  (and  the 
electric  field)  at  every  cell  (volume)  in  the  mesh.  The  BEM  requires  only  a  surface  mesh 
(although  volume  information  is  required  in  defining  the  different  “domains”)  and  calculates  the 
electric  potential  (and  the  electric  field)  only  at  boundary  and  interface  faces.  The  two  methods 
will  be  discussed,  followed  by  a  brief  discussion  comparing  the  methods.  The  benchmarks 
comparing  the  two  methods  are  provided  in  Chapter  5. 

3.1  Boundary  Element  Method  (BEM) 

BEMs  deal  with  “boundary  elements”  which  separate  “domains”.  The  boundary  elements  are 
defined  by  a  surface  grid  S.  A  domain  is  a  region  in  space  with  a  unique  £r  or  where  no  solution 
is  being  computed  (an  empty  domain).  BEM  can  solve  both  bounded  and  unbounded  problems. 
For  unbounded  problems  the  potential  at  infinity  is  assumed  zero  (i.e.,  (t>(°°)  —  0). 

Significant  enhancements  were  made  to  the  CFDRC  BEM  solver  to  obtain  O(N)  computation 
times  for  both  the  Laplacian  and  the  Navier  equations.  The  first  phase  of  development  involved 
the  implementation  of  a  stable  fast  boundary  element  solver/  library  for  Laplace/Poisson,  Stokes, 
and  Navier  equations  using  low  order  interpolation  and  multipole-like  expansions. 

The  BEM  solver  is  capable  of  solving  different  equation  types  including  Laplace/  Poisson, 
Navier,  Stokes,  and  Hemholtz  at  low  frequency.  Boundary  conditions  BEM  supports  include 
Dirichlet,  Neumann,  and  third  kind  (convective  for  Laplace  and  spring  support  for  Navier). 
Multi-domains,  each  with  different  set  of  properties,  are  supported.  Continuity  of  the  primary 
variable  (temperature,  electric  potential,  displacement)  as  well  as  the  secondary  variable  (fluxes 
and  tractions)  are  automatically  handled  across  the  interface.  A  new  method  (developed  within 
CFDRC)  is  used  to  evaluate  the  boundary  integrals.  The  method  is  similar  to  the  multipole 
method,  but  is  applicable  to  any  boundary  integral  equation,  including  the  3-D  Navier  equation. 
Sources  terms  (terms  typically  placed  on  the  right-hand  side  of  the  equations  and  are  not  a 
function  of  the  primary  variable  solved)  are  approximated  by  particles.  Fast  summations  of  the 
sources  terms  are  used.  The  time  complexity  is  0(log  M)  for  each  summation  of  the  contribution 
of  M  particles.  A  “matrixless”  generalized  minimal  residual  (GMRES)  method  is  used  to  solve 
the  linear  equation  systems  iteratively.  A  diagonal  preconditioner  is  used  everywhere  except  on 
the  interfaces,  where  a  2x2  submatrix  preconditioner  is  used.  The  overall  time  complexity  is 
O(N)  without  sources  and  0(N  log  M)  with  sources  for  each  iteration.  The  total  time  also 
depends  on  the  number  of  iterations  of  the  GMRES  method.  This  number  is  in  turn  problem 
dependent  and  preconditioner  dependent.  In  general,  about  N005  iterations  are  needed. 
Therefore,  the  time  dependence  is  roughly  O(N1,05~1,25)  or  0(N  '  log  M).  The  spatial 
dependence  is  0(N+M).  For  small  problems  (with  up  to  a  few  hundred  elements),  conventional 
0(N2)  BEM  can  be  faster  than  accelerated  BEM.  The  conventional  BEM  solver  can  be  activated 
by  a  simple  parameter  setting.  The  BEM  libraries  are  accessible  through  C  function  calls.  The 
user  can  print  the  results  in  CFD-VIEW  format  or  VRML. 
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3.1.1  Overview 

The  BEM  for  the  electric  problem  uses  the  “FastBEM”  solver.  The  FastBEM  solver  is  a  general 
purpose  BEM  solver  developed  at  CFDRC  and  capable  of  solving  Navier-Stokes  and  Laplace¬ 
like  equations.  In  general,  FastBEM  can  be  used  to  solve  any  equations  of  the  form 

J [H (r ' - r ) u (r ') - G (r ' -r) q (r ') ] dSr '  =  J[K(r'-r)s(r')dV(r'),r  e  f2  (3.1) 

s  V 

where  H,  G,  and  K  are  kernel  functions,  s  is  a  source  term,  u  is  the  primary  variable,  q  is  the 
“flux”  variable,  and  the  integrals  are  over  a  closed  surface  S  containing  a  volume  V.  The 
electrostatic  module  will  use  the  Laplace-like  form  of  Eq.  3.1  which  is 
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y4  nxjj 
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(3.2) 


where  r  =  |r'  -  rj,  n  is  the  normal  to  the  surface  S,  §  is  the  electric  potential,  8  is  the  electric 
permittivity,  and  p  is  the  space  charge. 

The  surface(s)  S  is(are)  divided  into  boundary  elements,  and  hierarchically  arranged.  For 
example,  the  elements  in  S  can  be  divided  into  subgroups  SI  and  S2  such  that  S  =  SI  +  S2. 
Similarly,  SI  =  Sll  +  S12,  S2  =  S21  +  S22,  and  so  on.  Any  tree  structure  can  be  used  for  the 
hierarchical  arrangement,  an  octree  is  used  in  many  multipole  methods.  FastBEM  uses  a 
balanced  axis  aligned  binary  tree. 

In  general  FastBEM  is  flexible  and  able  to  solve  different  equation  types.  The  penalty  for  this 
flexibility  is  more  memory  and  CPU  time.  Therefore,  a  technique  called  code  specialization 
was  introduced.  A  specialized  version  of  FastBEM  has  been  created  to  handle  only  a  single  type 
of  equation,  with  corresponding  gains  in  memory  and  time  efficiency  (by  more  than  100%).  If 
needed,  the  original  general  version  can  be  also  compiled,  with  a  simple  compiling  option. 

3.1.2  Space  Charge  Approximation 

The  space  charge  source  term  is  approximated  as  a  summation  of  point  charges  in  FastBEM. 
That  is, 


n  sources 

p(r)«  ZQi5(r-ri) 

i=l 


(3.3) 


Qi=  fp(r)dvi 

Vi 


(3.4) 


where  rj  and  Qi  are  the  location  and  the  strength  of  the  point  charge  i  respectively,  ^sources  is  the 
number  of  point  charges,  5  is  the  direct  delta  function,  and  v,  is  the  volume  over  which  a  point 
charge  i  has  been  approximated. 
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3.1.3  Kernel  Expansion  and  Fast  Integration 

When  the  minimum  distance  between  a  position  r  and  a  surface  S  is  large  compared  to  the 
distance  between  the  furthest  two  points  on  a  surface  S  the  kernel  functions  become  smooth  and 
can  be  approximated  by  low  order  expansion.  That  is  for  a  scaling  parameter  0  <  0  <  1.  In  the 
extreme  case  0  «  1, 

When  the  minimum  distance  between  a  position  r  and  a  surface  S  is  large  compared  to  the 
distance  between  the  furthest  two  points  on  a  surface  S  the  kernel  functions  in  8.1  become 
smooth  and  can  be  approximated  by  low  order  expansion.  That  is  for  a  scaling  parameter  0  <  0  < 
1. 


min  ]|r  —  r'|| 
r'eS 


max  ||r"-r'||| 
r',  r"  e  S 


(3.5) 


In  the  extreme  case  0  «  1,  a  kernel  can  be  approximated  by  a  single  term  approximation.  For  a 
kernel  K(r '  -  r)  the  approximation  is  a  single  value  of  K  at  the  center  rcenter  of  a  surface  S,  that  is 
K(r'  -  r)  «  K(rCenter  -  r).  For  0  =  0  no  (r,S)  pair  can  satisfy  Eq.  (8.5),  so  all  calculations  have  to 
be  performed  with  no  optimization.  In  this  case,  the  conventional  BEM  is  recovered.  Currently  m 
FastBEM  the  kernels  are  expanded  using  the  Taylor  series  expansion  if  they  meet  the  criteria  for 
expansion. 

For  quadpole  expansions  in  their  optimal  performance  range,  the  solution  error  m  primary 
variables  (temperature,  electric  potential  and  displacement)  is  typically  5%  to  10%.  The  error  in 
secondary  variables  (heat  flux,  electric  field,  and  traction)  is  even  larger.  Higher  order 
expansions  are  available.  The  expansions  have  been  implemented  in  a  hierarchy  of  C++  classes. 
Two  general  expansion  schemes  of  any  user  specified  order  are  available. 

Once  a  position  r  and  a  surface  S  in  Eq.  (3.2)  satisfies  the  distance  condition  in  Eq.  (3.5),  the 
kernel  can  be  approximated  by  a  few  terms  of  its  expansion  (such  as  Taylor  series),  and  the 
integral  can  be  quickly  evaluated.  If  a  position  r  relative  to  a  surface  S  does  not  satisfy  the 
distance  condition  in  Eq.  (3.5)  one  moves  down  the  sort  “tree”  to  subsets  of  the  surface  S  (SI, 
S2,  etc.)  to  see  if  (r,Sl),  (r,S2),  etc.  pairs  satisfy  the  distance  condition.  If  any  pairs  satisfy  the 
distance  condition  then  the  integral  is  approximated  by  a  set  amount  of  terms  of  a  senes 
expansion  (such  as  a  Taylor  series  expansion),  otherwise  the  surface  is  broken  down  further  and 

the  condition  is  tested  for  again. 

Even  more  efficiency  is  gained  when  r  is  also  arranged  into  hierarchies.  So,  the  criteria  for  the 
expansion  of  a  kernel  is  not  only  between  a  position  in  space  r  and  a  surface  S  (r,S)  but  (S,T) 
pairs  between  two  surface  groups  one  in  r  e  T  and  the  other  in  r'  €  S.  Whether  (T,S)  =  (T1,S)  + 
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(T2,S)  or  (T,S)  =  (T,S1)  +  (T,S2)  is  used  depends  on  which  grouping  increases  the  efficiency  of 
the  calculation. 

For  piecewise  constant  interpolations,  semi-analytical  integration  methods  are  usually  used  for 
all  kernels.  For  higher  order  interpolations,  pure  numerical  integration  was  used.  The  latter  is 
very  slow,  and  consumes  considerable  CPU  time.  During  testing  it  was  found  that  for  some 
problems  in  elastostatics  at  least  linear  interpolation  is  needed.  Therefore,  a  similar 
semi-analytical  integration  method  has  been  implemented  for  the  higher  order  interpolations 
(linear  to  cubic)  for  quadrilateral  and  triangle  boundary  elements.  For  other  polygon  elements, 
piecewise  constant  interpolation  is  retained.  Semi-analytical  integration  schemes  were  introduced 
in  the  BEM  model  for  the  singular  integrals  and  integrals  in  which  the  source  point  and  the 
integration  domain  are  co-planar.  This  capability  is  now  available  for  the  Laplacian/Poisson  and 
Navier  kernels  with  constant  and  linear  interpolations.  The  analytical  integrals  reduce  the 
computational  time  significantly,  e.g.,  for  a  capacitance  extraction  problem  (for  a  unit  cube  with 
2400  faces),  the  analytical  integrals  reduce  computation  time  by  a  factor  of  3,  i.e.,  2.87  s  from 
9.08  s  with  numerical  integration. 

3.1.4  FastBEM  Implementation 

FastBEM  is  written  in  C++  and  is  accessed  by  ACE+  through  C  subroutine  calls.  FastBEM  deals 
with  “boundary  elements”  which  separate  “domains”.  A  domain  is  a  region  in  space  with  a 
unique  £r  or  where  no  solution  is  being  computed  (an  empty  domain).  FastBEM  requires  a  list  of 
nodes  which  make  up  the  boundary  elements,  the  nodes  physical  position,  the  node  to  boundary 
element  connectivity,  a  list  of  domains,  the  domain  present  on  either  side  of  a  boundary  element, 
the  boundary  condition  on  a  side  opposite  the  end  of  an  empty  domain  if  present,  and  the 
location  and  strength  of  the  point  source  terms.  The  node  order  which  specifies  a  boundary 
element  follows  the  right-hand  rule  for  face  normal  definition.  A  boundary  element  face  “side” 
uses  the  normal  from  the  face  to  distinguish  between  sides. 

The  problem  specification  (grid,  boundary  conditions,  material  properties,  space  charge  volume 
conditions,  and  solution  method)  is  inputted  into  the  electrostatic  module  using  a  DTF  file 
format.  If  the  solution  method  is  the  BEM  then  the  problem  specification  data  is  reorganized  for 
submittal  to  FastBEM  subroutines. 

Rules  were  established  for  “faces”  which  will  not  be  submitted  to  FastBEM  as  boundary 
elements:  (i)  faces  between  two  non  -1  domains  with  a  boundary  condition  other  than 
dielectric/dielectric,  (ii)  faces  between  two  sr  =  0  (i.e.  “empty”,  or  metal)  domains,  (iii)  faces 
between  two  similar  domains  (sr  and  mat_type  the  same),  and  (iv)  faces  with  “ignore”  boundary 
condition. 

The  FastBEM  libraries  for  solution  of  Laplace-like  equations  are  used  within  the  CFD-ACE+ 
environment.  A  connection  of  the  FastBEM  Laplace-like  equation  solver  to  the  GUI  was 
implemented  through  C  subroutine  interface  calls  supplied  by  the  FastBEM  libraries.  Problem 
specification  data  such  as  boundary  conditions,  domain  information,  and  permittivity  values,  and 
charge  distributions  are  processed  into  a  form  recognized  by  FastBEM.  A  special  permittivity 
value  is  used  to  flag  to  FastBEM  the  presence  of  a  closed  domain. 
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3.1.5  FastBEM  Grid  Adaptation 

Grid  adaptation  capability  (for  BEM)  has  been  added  as  an  independent  module.  It  reads  in  the 
original  grid  file  and  the  initial  solution  by  FastBEM,  and  generates  a  new  grid  adapted  to  the 
solution.  Figure  3.1  demonstrates  grid  adaptation  in  a  BEM  solution. 


Figure  3-1.  Grid  Adaptation  in  FastBEM 


3.2  Finite  Volume  Method  (FVM) 

The  FVM  for  the  electric  problem  uses  finite  differencing  and  the  Divergence  Theorem  to 
discretize  the  solution  space 

jV.fdV  =  f(f*n)dS  (3.6) 

v  s 

cellfaces 

ZAjfj 

V.fcell-  “y  (3-7) 

The  normal  diffusive  fluxes  (f  •  n )  are  approximated  as  the  ratio  of  the  difference  of  cell  center 
values  and  cell  to  face  distances.  The  resulting  difference  equation  is  solved  using  an  iterative 
conjugate  gradient  solver  (CGS). 

3.3  BEM  vs.  FVM 

The  cpu  “price”  per  element  (cell)  is  much  greater  for  BEM  than  for  FVM.  A  typical  question 
then  is,  “when  is  one  method  “better”  or  “worse”  to  use  than  the  other?”  In  most  cases,  FVM  is 
the  “better”  method.  There  are  two  instances  when  BEM  is  the  better  or  only  choice.  One  is  if  the 
“volume  cell”  to  “boundary  +  interface  face”  ratio  is  on  the  order  of  100  or  1000.  In  this  case  the 
higher  cpu  cost  per  element  of  BEM  is  offset  by  the  fewer  number  of  elements  in  the  BEM 
computation  (boundary  +  interface  faces)  compared  to  the  number  FVM  will  use  (volume  cells). 
Another  instance  when  BEM  may  be  chosen  over  FVM  is  for  unbounded  problems,  that  is 
problems  with  an  infinite  domain.  The  FVM  requires  extra  cells  to  extend  the  boundary  to  a 
region  where  the  solution  no  longer  varies,  approximating  the  infinity  point.  With  BEM  no  extra 
meshing  is  required  to  extend  the  domain. 
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4. 


INTEGRATION  OF  PHYSICAL  MODELS  INTO  COMMERCIAL  CAD  TOOL 


4.1  CFD-ACE+MEMS  CAD  System 

The  main  objective  of  the  proposed  study  is  to  develop  an  integrated  commercial  CAD  tool  for 
MEMS  applications.  It  is  therefore  important  to  link  the  physical  models  with  other  peripheral 
software  (such  as  grid  generators,  geometry  modelers,  visualization  software,  Graphical  User 
Interfaces,  interfaces  to  other  CAD  packages,  etc.)  to  enhance  the  value  of  the  physical  models  as 
a  design  tool.  The  thermo-fluidic-mechanical-electric-magnetic  model  developed  during  this 
project  has  been  incorporated  into  a  design  software  environment,  the  CFD-ACE+MEMS  CAD 
system  shown  in  Figure  4-1.  The  design  environment  for  MEMS  which  has  been  developed 
includes  grid  generation  (CFD-GEOM),  data  visualization  (CFD-VIEW),  graphical  problem 
setup  (CFD-GUI),  and  coupled  fluidic,  thermal,  mechanical,  electric,  and  magnetic  physical 
model  solvers  (CFD-ACE,  CFD-ACE(U)). 


CFD-ACE  + 


Figure  4-1.  The  CFD-ACE+MEMS  CAD  System 

Geometry  Modeling  and  Grid  Generation  Tool:  The  CFD-GEOM  software  developed  by 
CFDRC  is  used  to  define  the  geometric  details  of  MEMS  devices.  CFD-GEOM  is  a  NURBS 
based  3-D  geometric  modeler.  It  has  points,  lines,  curves  and  surfaces  as  the  basic  geometric 
entities  from  which  complex  geometries  can  be  defined.  It  is  possible  to  import  geometric 
definitions  from  other  CAD  tools  (Pro/ENGINEER,  MICROSTATION,  ICEM,  and  AutoCAD) 
through  IGES  files.  GEOM  can  also  import  ECAD  formats  such  as  CIF  and  GDSII. 
Additionally,  GEOM  has  interfaces  to  facilitate  geometry  and  data  transfer  from  other 
commercial  CAD  packages  such  as  Pro/ENGINEER,  AutoCAD. 

CFD-GEOM  can  be  used  to  create  both  structured  and  unstructured  meshes  for  the  simulations. 
The  basic  entities  for  creating  meshes  include  edges,  faces,  loops,  surfaces,  blocks,  and  domains. 
It  is  also  possible  to  form  prisms  (triangles  extruded  in  the  third  dimension)  which  are  useful  for 
device  geometries  which  have  most  of  the  spatial  variation  in  two  dimension. 

Visualization  Tool:  CFD-VIEW,  developed  by  CFDRC,  is  used  to  visualize  the  results  of 
simulations.  The  input  formats  recognized  by  this  package  include  FAST-U,  hybrid,  mixed-u, 
polyhedral,  PLOT3-D  data  sets  (NASA  Standard)  or  the  DTF  format  developed  by  CFDRC. 

Visual  Environment:  The  Visual  Computing  Environment  (V CE)  capability  being  developed  at 
CFDRC  for  parallel  execution  of  multi-disciplinary  analysis  codes  will  also  be  available  in  the 
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CAD  tool  environment.  VCE  is  an  object-oriented  environment  for  running  multiple  software 
modules  on  a  network  of  heterogeneous  workstations. 

Graphical  User  Interfaces  (GUIs):  Complete  access  to  the  CAD  tool  is  provided  through  a 
master  GUI.  This  GUI  facilitates  access  to  any  of  the  individual  modules  such  as  CFD-ACE, 
CFD-GEOM,  CFD-VIEW,  etc. 

Physics  Solvers:  CFD-ACE  and  CFD-ACE(U)  are  modules  of  the  CFD-ACE+MEMS  analysis 
package.  They  provide  the  tools  to  specify,  and  iteratively  solve,  the  equations  that  describe  the 
physical  processes  for  fluidic,  thermal,  mechanical,  electric,  and  magnetic  problems.  The 
physical  models  are  solved  on  a  2-D,  cylindrical,  or  3-D  multi-domain  grid.  The  grid  can  be 
structured,  unstructured,  hybrid,  or  adaptive  Cartesian.  There  is  also  a  moving  grid  capability  to 
track  moving  and  deforming  bodies  and  surfaces.  Strong  links  exist  with  data  transfer  facilities 
(DTF),  parallel  execution  on  multiple  machines  using  the  message  passing  interface  (MPI),  and  a 
link  with  ADAPTER  for  solution  based  grid  refinement/adaptation.  The  physical  models  are 
implemented  in  a  highly  modularized  code  architecture  which  facilitates  the  addition  of  future 
physical  models. 

In  addition  to  the  above  mentioned  existing  software  modules,  a  virtual  component  library  and  a 
MEMS  knowledge  base  (to  be  developed  under  an  ongoing  project)  will  be  integrated  with  the 
software  environment.  The  component  library  will  contain  a  set  of  predefined,  parameterized 
components  which  can  be  used  to  design  a  complicated  system.  The  knowledge  base  will  contain 
flow  and  process  parameters  and  their  relationships  for  these  components.  It  will  also  contain 
design  rules  and  guidelines  for  integration  with  other  components. 

4.2  Electric  Model  in  CFD-ACE+MEMS  CAD  System 

The  electric  model  physics  solver  was  discussed  in  the  previous  chapters.  In  this  section  the 
electric  problem  specification  capabilities  in  CFD-ACE+MEMS  will  be  discussed.  First,  the 
problem  specifications  that  are  common  to  both  FVM  and  BEM  solution  methods  of  the  electric 
model  are  given.  The  common  problem  specification  options  are  the  boundary  conditions, 
relative  permittivity  material  property,  and  the  space  charge  volume  condition.  This  is  followed 
by  a  discussion  of  additional  problem  set-up  specific  to  FVM  and  BEM  solution  methods. 
Finally,  the  variables  deposited  to  the  DTF  formatted  input/output  file  will  be  discussed. 

A  user  can  choose  between  FVM  and  BEM  using  the  dialog  shown  in  Figure  4-2. 
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Options 


J  Incompressible  Flow 
J  Compressible  Flow 
J  Turbulence 
J  Heat  Transfer 
J  Radiation 
j  Mixing 
J  Reaction 
J  Deformation 
J  Stress 
■  Electrostatics 
J  BEM  Solver 
'•  FVM  Solver 

J  Free  Surfaces  (VOF) 
-j  Curing  Chemistry 
-i  Surface  Tension 
j  Voids 

J  Plasma 

J  Electromagnetics 


Figure  4-2.  Dialog  Where  Electric  Model  is  Chosen  and  Either  FVM  or  BEM  Solver  is  Selected 


4.2.1  Electrostatic  Boundary  Conditions 

The  five  boundary  conditions  available  in  the  module  and  shown  in  Figure  4-3  are  specifying  the 
value  of  the  potential  at  a  cell  face  (Eq.  (4.1)),  specifying  the  negative  normal  potential  gradient 
(electric  field  normal)  at  a  cell  face  (Eq.  (4.2)),  specifying  a  mixture  of  fixed  potential  and 
electric  field  normal  (Eq.  (4.3)),  indicating  a  cell  face  is  at  a  dielectric  interface  which  is  useful 
when  specifying  an  unbounded  problem,  and  an  ignore  boundary  condition  which  is  useful  when 
multi-disciplinary  problems  are  specified. 
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Figure  4-3.  Electrostatic  Boundary  Conditions  Window  at 
Model. “Bond.Cond.”.“Surface  BC”. Values 

4.2.1.1  Fixed  Potential  or  “Potential  tVoItageV’.  The  “Potential  (Voltage)”  boundary 
condition  in  the  GUI  allows  the  specification  of  a  fixed  potential  at  a  cell  face.  Clicking  on 
“Potential  (Voltage)”  causes  a  “Voltage:”  labeled  box  to  appear.  One  then  specifies  the  voltage 
in  units  of  volts  that  the  selected  face(s)  should  be  fixed  at. 

4>  =  <f)0  =  Voltage  (4.1) 

4.2.1.2  Fixed  Flux  or  “Electric  Field  Normal”.  The  “Electric  Field  Normal”  boundary 
condition  in  the  GUI  allows  the  specification  of  a  fixed  electric  field  normal  at  a  cell  face 
weighted  by  a  relative  permittivity.  Clicking  on  “Electric  Field  Normal”  causes  a  “Flux:”  labeled 
box  to  appear.  One  then  specifies  the  electric  field  normal  or  flux  in  units  of  volts/meter  that  the 
selected  face(s)  should  be  fixed  at. 


-srV<|)  ■n  =  Flux  (4.2) 

The  normal  direction  is  away  from  the  mesh  at  mesh  end,  toward  the  solid  at  a  solid/fluid 
interface,  and  dependent  on  the  node  order  of  the  face  at  a  solid/solid  interface. 

4.2.1.3  Mixture  of  Fixed  Potential  and  Flux  or  “Mixed”.  The  “Mixed”  boundary  condition  in 
the  GUI  allows  the  specification  of  a  mixture  of  fixed  potential  and  fixed  electric  field  normal  at 
a  cell  face.  Clicking  on  “Mixed”  causes  three  boxes  labeled  “a:”,  “b:”,  and  “c:”  to  appear.  One 
then  specifies  the  values  for  the  constants  a,  b,  and  c  which  fit  the  following  equation. 


a<j>-berV<j)-h=c  (4.3) 

4.2. 1.4  Dielectric/Dielectric  Interface  or  “Dielectric”.  The  “Dielectric”  boundary  condition  in 
the  GUI  is  used  to  indicate  that  a  dielectric/dielectric  interface  exists.  Interfaces  which  separate 
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regions  of  different  permittivity  within  the  grid  are  dealt  with  automatically  in  the  code  so  there 
is  no  need  to  apply  the,  “Dielectric”  boundary  condition  at  these  interfaces.  The  boundary 
condition  is  used  for  electrostatic  problems  when  the  BEM  solver  is  used.  The  boundary 
condition  allows  the  user  to  include  faces  in  the  calculation  where  the  voltage  value  is  desired 
but  no  knowledge  of  the  potential  or  its  gradient  is  known.  Clicking  on  “Dielectric”  causes  no 
labeled  boxes  to  appear.  This  boundary  condition  does  not  require  any  constants  to  be  specified. 

4.2.1.5  No  Electrostatic  Boundary  Condition  Exists  or  “Ignore”.  The  “Ignore”  boundary 
condition  in  the  GUI  is  used  to  indicate  that  a  set  of  faces  has  no  electrostatic  boundary 
conditions.  This  boundary  condition  is  useful  when  multi-disciplinary  problems  are  being 
specified.  A  set  of  faces  may  require  a  boundary  condition  to  be  specified  in  one  problem  type 
but  have  no  boundary  conditions  applicable  for  another  problem  type.  The  “Ignore”  boundary 
condition  allows  one  to  specify  that  a  set  of  faces  has  no  electrostatic  boundary  condition 
although  it  may  have  a  boundary  condition  for  another  problem  type.  Clicking  on  “Ignore” 
causes  no  labeled  boxes  to  appear.  This  boundary  condition  does  not  require  any  constants  to  be 
specified. 

4.2. 1.6  Symmetry.  The  symmetry  boundary  condition  in  the  GUI  is  used  to  specify  areas  of 
symmetry  in  the  grid.  That  is,  specify  symmetry  planes.  Mathematically  it  is  equivalent  to 
-srV(|)-h=0  or  the  field  normal  at  a  face  is  equal  to  zero.  This  boundary  condition  does  not 
require  any  constants  to  be  specified.  It  should  be  noted  that  for  unbounded  problems  solved  with 
BEM  there  is  no  possibility  of  specifying  a  symmetry  plane.  An  infinite  plane  would  need  to  be 
specified  with  a  zero  normal  field. 

4.2.2  Electrostatic  Volume  Conditions 

The  two  electrostatic  volume  conditions  that  can  be  specified  are  the  relative  electric  permittivity 
sr,  and  space  charge  p,  as  defined  in  Eq.  (2.1). 


Solid  Properties 


-Materials  - 


I-  Electrostatics  - 


Relative  Permittivity  111 


Space  Charge  ij0 


Close] 


Figure  4-4.  Electrostatic  Volume  Conditions  Window  at  Model.Prop.Solid.Electrostatics  or 
Model.Prop.Gas.Electrostatics 
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4.2.2.1  Electric  Permittivity.  Only  the  relative  permittivity  £r  of  the  total  electric  permittivity  . 
sr£0  is  set  in  the  GUI.  The  value  used  for  the  permittivity  of  free  space  is  £0  =  8.854187  xlO12 
F/m  or  C2/Nm2.  The  relative  permittivity  defaults  to  1.0  for  all  cells  in  the  grid. 

To  specify  a  metal  (or  an  empty  domain  in  FastBEM)  the  relative  electric  permittivity  for  the 
group  of  cells  should  be  set  to  zero.  No  electric  field  exists  in  a  metal  and  a  relative  electric 
permittivity  of  zero  is  used  in  the  code  as  a  flag  to  indicate  that  a  group  of  cells  should  not  be 
included  in  the  solution  to  Poisson’s  equation.  In  the  FVM  the  cells  are  “blocked”  so  that  the 
electric  potential  and  field  are  not  solved  for  those  cells.  In  the  BEM,  setting  a  group  of  cells  to 
zero  relative  electric  permittivity  is  equivalent  to  establishing  an  empty  domain.  If  fixed 
potential,  fixed  flux,  or  mixed  boundary  conditions  are  set  on  a  face  one  of  the  two  domains  the 
face  separates  must  be  an  empty  domain. 

4.2.2.2  Space  Charge.  The  space  charge  p  of  a  group  of  cells  can  be  specified  in  the  GUI  in 
units  of  C/m5.  The  space  charge  defaults  to  0.0  for  all  cells  in  the  grid. 

4.2.3  FVM  Control  Parameters 

The  FVM  control  parameters  for  the  electrostatic  problem  are  accessed  at  different  places  in  the 
GUI.  The  iteration  number  and  final  residual  are  accessed  at  Solve.Control. Electrostatics.  The 
solver,  relaxation,  and  limits  are  accessed  at  Solve.Control.  The  FVM  control  parameters 
windows  are  shown  in  Figure  4-5. 
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Figure  4-5.  An  Overview  of  the  FVM  Controls.  The  control  menu  (upper  left)  is  used  to  reach 
the  solver,  relaxation,  limits,  iteration,  and  residue  controls. 

4.2.3.1  Maximum  Iteration  Number 

>0 

The  “Maximum  Iteration  Number”  controls  the  maximum  number  of  iterations  to  be  made  by  the 
FastBEM  algorithm  before  exiting.  The  default  value  is  200. 

4.23.2  Residue 


>0 

The  “Residue”  is  the  desired  residue  reduction  factor  for  the  electric  potential.  The  residue  is 
normalized  to  one  at  the  beginning  of  the  calculation.  The  calculation  will  stop  when  either  the 
desired  electric  potential  residue  has  been  achieved  or  the  maximum  iterations  have  been 
exceeded.  The  default  value  is  1.0  xlO'4. 
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4.2.3.3  Electric  Potential  Limits.  A  maximum  and  minimum  value  can  be  set  on  the  value  of 
the  electric  potential  as  it  is  being  solved. 

4.2.3.4  Under  Relaxation.  Under  relaxation  can  be  used  when  solving  the  electric  potential 
using  FVM  which  weights  the  current  solved  value  with  the  previous  solution. 

4.2.3.4  Linear  Solver.  The  type  of  linear  solver  can  be  chosen  (conjugate  gradient,  arithmetic 
multi-grid),  as  well  as  the  maximum  number  of  sweeps  to  take  in  the  solver  and  the  residual 
reduction  criteria  for  exiting  the  solver. 

4.2.4  BEM  Control  Parameters 

The  BEM  control  parameters  for  the  electrostatic  problem  are  accessed  at 
Solve.Control.Electrostatics.  The  Electrostatic  option  is  not  available  in  the  Control. Solve  menu 
unless  the  Model.Problem_Type.Options. Electrostatics. BEM  Solver  option  has  been  selected. 
The  BEM  control  parameters  window  is  shown  in  Figure  4-6. 


Electrostatics:  BEM  Solver 


Interpolation  Order  |[0 
Maximum  Iteration  Number  H200 
Gmres  Restart  Steps  jllO 
Expansion  Mode  j[l 
Expansion  Order  |[2 
Distribution  Criteria  ijO .  5 

Estat  Residue  |]0 . 000 1 
J  Unbounded 


Default 


Electrostatics:  BEM  Solver 


Interpolation  Order 

,P 

Maximum  Iteration  Number 

.[200 

Gmres  Restart  Steps 

!io 

Expansion  Mode 

H 

Expansion  Order 

[2 

Distribution  Criteria 

j0. 5 

Estat  Residue 

[0.0001 

jv  Unbounded)  Er 

[1 

OK  I  Default 


Figure  4-6.  BEM  Control  Parameters  Window  at  Solve.Control.Electrostatics 

4.2.4.1  Interpolation  Order 

0  -  uniform  (anything  larger  than  0  is  very  computational  expensive) 

1  -  linear 

The  “Interpolation  Order”  controls  the  interpolation  level  on  a  face  where  the  potential  and/or 
gradient  will  be  solved.  Higher  interpolation  order  leads  to  better  accuracy  in  results  but 
increases  the  computational  expense  of  the  calculation  drastically.  Only  zeroth  (uniform)  and 
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first  order  (linear)  interpolation  are  currently  available.  In  addition  first  order  interpolation  is 
only  available  for  n-polygon  faces  of  n  =  3  or  4  (i.e.,  triangles  and  quadrilaterals)  faces.  If  first 
order  interpolation  is  specified  any  n-sided  polygon  faces  greater  than  four  will  use  zeroth  order 
interpolation.  Zeroth  order  introduces  only  one  node/face  at  face  center.  First  order  interpolation 
introduces  n  nodes  for  an  n-sided  polygon  face.  Inputted  values  less  than  zero  are  changed  to 
zero.  Inputted  values  greater  than  1  are  changed  to  1 .  The  default  value  is  0. 

4.2.4.2  Maximum  Iteration  Number 


>0 

The  “Maximum  Iteration  Number”  controls  the  maximum  number  of  iterations  to  be  made  by  the 
FastBEM  algorithm  before  exiting.  The  default  value  is  200. 

4.2.4.3  Gmres  Restart  Steps 

>0 

A  matrix-less  GMRES  method  is  used  to  solve  the  equation  systems  iteratively.  A  diagonal 
preconditioner  is  used  everywhere  except  on  the  interfaces,  where  a  2x2  submatrix 
preconditioner  is  used.  The  “Gmres  Restart  Steps”  parameter  controls  the  number  of  sweeps 
within  the  linear  solver  for  each  iteration.  The  default  value  is  10. 

4.2.4.4  Expansion  Mode 
1 

The  “Expansion  Mode”  controls  the  method  of  approximating  the  kernel  function  as  a  series 
expansion.  Currently  only  Taylor  series  expansion  of  the  kernel  function  is  available.  Multipole 
expansion  itself  is  not  implemented  due  to  its  lake  of  generality.  In  the  future  other  expansion 
modes  may  become  available.  The  default  and  only  possible  value  is  1  so  any  changes  made  here 
will  not  effect  the  solution. 

4.2.4.5  Expansion  Order 


0-4 


The  “Expansion  Order”  is  the  order  of  the  kernel  series  expansions  (currently  only  the  Taylor 
series  expansion  is  available).  The  default  value  is  2. 

4.2.4.6  Distribution  Criteria 


0.0<x<  1.0 

The  “Distribution  Criteria”  is  the  multipole  accuracy  parameter.  Loosely  put  the  parameter 
controls  the  level  of  “fast”  in  the  FastBEM  algorithm  that  is  used.  The  “Distribution  Criteria”  is 
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connected  to  the  0  of  Eq.  (3.3).  A  value  of  zero  is  conventional  BEM  with  no  attempts  made  to 
speed  up  the  algorithm.  The  solution  will  have  a  high  level  of  accuracy  but  the  solution  speed 
will  be  very  slow.  As  the  parameter  is  increased  from  zero  the  level  of  optimization  increases  as 
attempts  are  made  to  speed  the  level  of  convergence.  Mechanisms  such  as  clumping  together 
boundary  elements  into  an  equivalent  boundary  element  are  used.  The  benefits  of  optimization 
(larger  Distribution  Criteria  values)  are  faster  solutions.  The  drawbacks  are  a  lower  level  of 
accuracy  of  the  solutions  and  an  increase  in  the  memory  requirements  to  achieve  a  solution.  The 
suggested  value  is  around  0.5.  Inputted  values  greater  than  one  or  less  than  zero  are  changed  to 
0.5.  The  default  value  is  0.5. 

4.2.4.7  BEM  Residue 


>0 

The  “BEM  Residue”  is  the  desired  residue  reduction  factor  for  the  electric  potential.  The  residue 
is  normalized  to  one  at  the  beginning  of  the  calculation.  The  calculation  will  stop  when  either  the 
desired  electric  potential  residue  has  been  achieved  or  the  maximum  iterations  has  been 
exceeded.  The  default  value  is  1.0  xl0‘4. 

4.2.4.8  Unbounded.  The  “Unbounded”  toggle  specifies  if  the  problem  is  bounded  or 
unbounded.  A  bounded  problem  assumes  the  boundary  elements  bound  the  solution  space.  An 
unbounded  problem  assumes  a  solution  is  unbounded  with  the  electric  potential  <J>(oo)  =  0. 
“Unbounded”  is  true  (toggle  on)  if  the  problem  is  an  unbounded  problem,  false  (toggle  off)  if  the 
problem  is  bounded. 

4.2.4.9  Relative  Permittivity  sr 

>  0.0 

The  relative  permittivity  in  the  unbounded  region.  The  relative  permittivity  in  the  unbounded 
region  can  not  be  specified  unless  the  unbounded  option  has  been  selected.  The  default  value  is 
1.0. 


4.2.5  Electric  Model  DTF  Output 

The  values  deposited  by  the  electrostatic  module  into  the  DTF  file  are  different  depending  on  the 
method  used  (BEM  or  FVM)  to  solve  Poisson’s  equation. 

4.2.5.1  BEM  Electrostatic  Output 


cfd-view 

description 

symbol 

units 

elpot 

electrostatic  potential 

<t> 

V 

epsr_e_En 

electric  field  normal 

-£rV<t)-h 

V/m 

pestat 

electrostatic  pressure 

P 

N/m2 
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All  three  values  are  calculated  only  on  boundary  or  interface  faces  and  averaged  onto  the  vertices 
using  1/r  weighting.  As  a  result  values  at  vertices  at  the  comer  of  the  grid  may  differ 
significantly  from  their  original  face  values.  The  electric  field  normal  weighted  by  the  relative 
permittivity  uses  the  normal  local  to  the  face.  The  electrostatic  pressure  force  is  the  normal  force 
directed  toward  an  empty  domain.  Thus  a  positive  pressure  pushes  the  face  “toward”  the  empty 
domain  and  a  negative  pressure  force  would  pull  a  face  “away”  from  the  empty  domain 

4.2.5.2  FVM  Electrostatic  Output 


cfd-view 

description 

symbol 

units 

elpot 

electrostatic  potential 

<t> 

V 

efieldx 

electric  field  x  component 

Ex 

V/m 

efieldy 

electric  field  y  component 

Ey 

V/m 

efieldz 

electric  field  z  component 

Ez 

V/m 

pestat 

electrostatic  pressure 

P 

N/m2 

All  values  are  calculated  at  cell  centers  and  averaged  onto  the  vertices  using  1/r  weighting.  As  a 
result  values  at  vertices  at  the  comer  of  the  grid  may  differ  significantly  from  their  original  cell 
values.  The  three  components  of  the  electric  field  E  =  Exx  +  Eyy  +  Ezz  are  given  as  the  three 
components  of  the  cartesian  coordinate  system.  At  material  interfaces  (faces  where  the  electric 
permittivity  changes)  the  electric  field  is  discontinuous.  This  is  not  seen  in  CFD-VIEW  where 
the  cell  center  values  have  been  averaged  onto  cell  vertices  before  output  to  the  DTF  file. 
Currently  there  is  no  way  to  assign  two  different  values  to  the  same  vertex.  The  electrostatic 
pressure  force  is  the  normal  force  directed  toward  an  empty  domain.  Thus  a  positive  pressure 
pushes  the  face  “toward”  the  empty  domain  and  a  negative  pressure  force  would  pull  a  face 
“away”  from  the  empty  domain 

4.3  Prototype  for  Web  Based  Interface 

An  HTML- VRML  based  interface  was  developed  to  allow  job  submission  to  FastBEM  via  the 
Web.  The  user  specifies  as  input  a  simple  CSG  model,  as  well  as  parameters  relating  to  mesh 
size,  problem  type,  and  boundary  conditions.  The  server  then  generates  the  grid,  solves  the 
problem,  and  returns  the  results  to  the  user's  browser,  along  with  a  link  to  the  3-D  plot  in  VRML 
format.  Sample  results  are  shown  in  Figure  4-7  to  Figure  4-10. 
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Figure  4-7.  Snapshot  of  Upper  Portion  of  GUI  for  CFD-FastBEM 
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Figure  4-8.  Snapshot  of  Lower  Portion  of  GUI  for  CFD-FastBEM 


Figure  4-9.  Results  from  On-line  Execution  of  FastBEM 
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Figure  4-10.  Visualization  of  FastBEM  Results  Using  VRML 
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5. 


VALIDATION  STUDIES  AND  STUDIES  OF  INDUSTRIAL  MICROSYSTEMS 


5.1  Benchmarks  of  BEM 

In  order  to  verify  the  0(N)  performance  of  the  developed  BEM  solver  several  benchmarking 
tests  were  carried  out  for  some  standard  problems  for  which  analytical  solutions  exist. 

The  first  problem  is  that  of  heat  conduction  inside  a  sphere  which  is  subjected  to  convection  and 
irradiation  on  the  surface.  The  temperature  distribution  is  shown  in  Figure  5-1.  The  second 
problem  is  that  of  a  cylindrical  vessel  deformation  under  surface  pressure.  The  displacement  is 
shown  in  Figure  5-2. 

Tables  5-1  and  5-2  summarize  the  time  and  memory  required  by  BEM  as  a  function  of  the 
number  of  elements,  N,  for  the  above  two  problems.  For  comparison,  the  performance  of  the 
conventional  BEM  (i.e.,  theta  =  0.0)  is  also  included. 

As  the  numbers  in  the  tables  show,  the  actual  performance  of  the  BEM  solver  is  indeed  very 
close  to  0(N)  for  each  iteration.  The  overall  time  performance  is  worse  for  the  second  problem 
(0(NA1.16)  to  0(NA1.32),  because  the  number  of  GMRES  iterations  grows  as  N  increases,  but 
this  should  be  expected  for  all  linear  solvers  of  the  CG  class. 


Figure  5-1.  Temperature  Distribution  on  a  Sphere  Subjected  to  Surface  Irradiation  and 
Convection 
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Figure  5-2.  Displacement  of  a  Cylinder  Under  Surface  Pressure 

Table  5-1.  Benchmark  Results  for  Heat  Conduction  Inside  a  Sphere 


Table  5-2.  Benchmark  Results  for  Deformation  of  a  Cylinder  Under  Surface  Pressure 


Additional  validation  of  the  FastBEM  model  was  carried  out  on  benchmark  problems  calculated 
by  NASTRAN’s  boundary  element  method  package,  TEAM.  Figure  5-3  shows  a  sample  list  of 
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3.7  Arc  Beam 


3.8  Plate  under  Gravity 


3.10  Spinning  Disk 


Figure  5-3.  Test  Cases  Used  by  TEAM 


Table  5-3.  Benchmark  Comparisons  with  TEAM 


_ _ _ 

Tapered  Plate  under  Edge  [  , 


Shear 


26.9  MPa 


Plate  with  a  Central  Crack  j  6.0066e-6  in 
j  Simple  Cantilever  Beam  0.33333e— 3  ini 


90°  Arc  Beam  Bending 

Tapered  Cantilever  under 
Gravity 

Flat  Spinning  Disk 


1.442e-3  in 


27.1  MPa 

4.13e-6  in 
-0.3510e-3  in 
1.768e-3  in 
1.584e-3  in 


":r~ 


-0.200  MPa  j  -0.200  MPa 
I.064e-3  in  i  1.066e-3  in 


5.2  Benchmarks  of  FVM  and  BEM  Applied  to  Electrostatics 


24.94  MPa  | 

198 

bilinear 

6.2704e-6  in  [ 

240 

bilinear 

-0.32535e-3  inf 

64 

bilinear 

1.927e-3  in  j 

48 

bilinear 

1.689e-3  in  j 

192 

-0.193  MPa  j 

70 

bilinear 

1.0Se-3  in  j 

208 

constant 

The  performance  of  FVM  and  BEM  were  compared  when  solving  Poisson’s  equation  on 
standard  test  problems.  Initial  testing  shows  FVM  faster  than  BEM  (for  simple  problems  such  as 


a  parallel  plate  capacitor)  and  more  accurate.  BEM  is  better  at  unbounded  problems  and 
problems  where  the  cells/boundaries  ratio  is  on  the  order  of  1000. 


5.2.1  Parallel  Plates 

The  first  benchmark  is  a  rectangle  with  two  different  electric  permittivity  materials  sandwiched 
between  square  parallel  contacts  as  shown  in  Figure  5-4.  It  is  a  bounded  problem  with  applied 
potentials  of  (J>(x  =  d)  =  <t>0  V,  <()  (x  =  0)  =  0  V.  All  other  surfaces  have  zero  normal  potential 
gradient  (i.e.,  -£rV<j)-h=0). 


Figure  5-4.  Parallel  Plate  Geometry 
The  electric  potential  at  x  =x0  is 


<KXo)  =  /  ,  £2Xf- -  <t>0 

s1(d-x0)  +  £2x0 


(5.1) 


The  uniform  electric  fields  at  either  side  of  the  x  =  x0  interface  are 


E,  = 


-e2 


1  £1(d-XG)  +  £2X( 


<t>0x 


(5.2) 


Eo  = 


£1(d-x0)  +  £2x0 


<f>o* 


(5.3) 


The  capacitance  is 


C  = 


A£i£2 


£l(d-X0)  +  £2X0 


(5.4) 


This  benchmark  was  done  on  a  DEC  500  with  d  =  10,  A  =  25,  and  <j)0  =  50  V  for  two  different 
cases  where  the  permittivity  of  the  two  regions  is  varied.  In  the  first  case  they  both  equal  the 
permittivity  of  free  space,  in  the  second  one  of  the  regions  has  a  permittivity  four  times  that  of 
free  space.  The  electric  potential  at  the  interface,  the  electric  field  in  each  region,  the 
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capacitance,  memory  requirements/  and  solution  time  were  compared  for  BEM  and  FVM.  As 
expected,  the  FVM  performs  much  better  than  BEM  for  this  case  in  accuracy,  speed,  and 
memory  requirements.  It  is  expected  that  BEM  will  be  more  approximate  given  that  accuracy  has 
been  sacrificed  to  gain  solution  speed.  The  BEM  predicted  capacitance  values  are  fairly  close  to 
theory  which  implies  that  on  average  the  calculated  field  values  at  the  contacts  is  accurate. 

(1)  Ei  =  e2  =  e0 


4)(x  =  5)  (V) 

Eix  (V/m) 

E2x  (V/m) 

C  (Farads) 

Mem 

CPU 

Theory 

25.00 

-5.00 

-5.00 

2.214x10'" 

MM:SS 

FVM 

25.00 

-5.00 

-5.00 

2.214x10'" 

3.5  M 

00:02 

BEM 

24.88  to 

25.17 

-5.25  to 
-4.85 

-5.25  to 
-4.72 

2.183x10'" 

18.0  M  . 

00:22 

(2)  Ei  =  e0,  e2  =  4e0 


<Kx  =  5)(V) 

E,x  (V/m) 

E2x  (V/m) 

C  (Farads) 

Mem 

CPU 

Theory 

40.00 

-8.00 

-2.00 

3.542x10'" 

MM:SS 

FVM 

40.00 

-8.00 

-2.00 

3.542x10'" 

3.6  M 

00:05 

BEM 

40.51 

-7.66  to 
-8.49 

-1.92  to 
+2.12 

3.332x10'" 

21.0  M 

00:23 

5.2.2  Concentric  Spheres 

Another  benchmark  case  is  two  spheres  with  different  radii  a  and  b  shown  in  Figure  5-5.  Electric 
potentials  of  <J)(r  =  a)  =  <|>0  V  and  <|)  (r  =  b)  =  0  V  are  applied  to  the  spheres.  The  potential  and 
electric  field  normal  were  calculated  at  the  center  of  25516  tetrahedrals  for  FVM  and  1128 
triangle  faces  for  BEM. 


Figure  5-5.  Concentric  Sphere  Geometry 
The  electric  potential  as  a  function  of  z  is 


<J)(z) =— <|)0  0<z<c 

c 


(5.5) 
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The  uniform  electric  fields  in  the  two  regions  are  equal  and  are  given  by 


Ei  =  E2 


The  capacitance  of  the  geometry  is  given  by 

Q  _  £  1 A 1  +£2A2 
c 


(5.6) 


(5-7) 


DEC  500,  a  =  1,  b  =  3,  <t>0  =  60  V,  e  =  8o 


<Kx  =  2)  (V) 

Er  (r=l)V/m 

Er  (r=3)V/m 

C  (Farads) 

Mem 

CPU 

Theory 

15.00 

90.00 

10.00 

1.669xlO"10 

MM:SS 

FVM 

15.00 

85.00 

11.00 

1.651xl0'10 

16  M 

00:15 

BEM 

n.a.  * 

85.0  to  90.1 

9.4  to  10.2 

1.648xlO'10 

51  M 

00:30 

*BEM  calculated  values  only  for  boundary  faces  at  r=l,3;  does  not  have  any  values  at  r=2 

The  values  in  the  table  are  approximate  due  to  the  amount  of  value  averaging.  First  there  is 
averaging  from  volume  center  to  nodes,  than  in  CFD-VIEW  where  the  point  and  line  probe 
average  from  the  node  values.  The  capacitance  calculation  which  is  done  internally  with  actual 
data  is  a  good  measure  of  how  close  the  average  electric  field  values  are  to  the  theory.  The  FVM 
field  values  are  off  because  the  cell  center  gradients  are  not  being  interpolated  properly  to  the 
faces.  The  values  reflect  the  value  at  cell  centers  near  the  boundary,  not  the  value  at  the 
boundary. 

5.2.3  Sphere  of  Space  Charge 

The  third  benchmark  study  of  FVM  and  BEM  is  the  electric  potential  and  field  due  to  a  sphere  of 
uniform  space  charge.  The  sphere  has  a  uniform  space  charge  of  p  and  a  radius  R  as  shown  in 
Figure  5.6.  The  total  charge  of  the  sphere  is  thus  Q  =  p(4/3)rcR3.  It  is  an  unbounded  problem  with 
the  potential  at  infinity  assumed  to  be  zero  (i.e.,  <|)(co)  =  0.0).  For  the  FVM  case  a  square  with 
length  300  was  generated  around  the  sphere  with  a  b.c.  of  <j)  =  0.0  V  at  the  six  square  sides  as  an 
artificial  truncation  point  for  the  domain.  In  theory  4>(r=300)  =  9/300  =  0.03  V.  Applying  b.c  of 
0.0  V  at  points  that  are  actually  0.03  V  introduces  error.  For  the  BEM  case  the  unbounded 
problem  needs  no  modification  of  the  mesh.  The  mesh  for  FVM  is  64996  tetrahedrals,  for  BEM 
the  potential  and  electric  field  normal  were  calculated  on  652  faces. 
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a 


Figure  5-6.  Sphere  of  Space  Charge  of  Radius  R  in  an  Optional  Box  with  Aides  of  Length  a  >  R 
The  electric  potential  as  a  function  of  the  radius  r  is 


m= 


4ne0r 


Q 


8tt80R 


,2  ^ 


3  - 


R" 


r>R 


r<R 


(5.8) 


The  electric  fields  in  and  out  of  the  sphere  as  a  function  of  r  are  given  by 


E(r)  = 


Q 

47ie0r2 


A 

r 


Qr  „ 
— — - r 

47te0R3 


r>R 


r<R 


(5.9) 


DEC  500,  R  =  1,  (a  =  300  for  fvm  case),  p  =  27 e0,  e  =  e0 


<t>(r=3)  (V) 

Er(r=3)  V/m 

Mem 

CPU 

Theory 

1.00 

MM:SS 

FVM 

2.9 

0.95  to  1.05 

40  M 

02:00 

BEM 

2.7  to  2.98 

0.81  to  1.05 

53  M 

00:45 

Neither  method  will  be  exact  given  that  the  sum  of  the  volumes  of  all  of  the  tetrahedrals  making 
up  the  meshed  sphere  will  not  equal  (4/3)tiR3.  The  FVM  accuracy  would  increase  if  the 
bounding  box  were  made  larger.  The  BEM  was  faster  than  FVM  and  as  accurate.  The  memory 
requirements  are  still  comparable.  The  set-up  time  for  the  BEM  was  also  much  less  than  FVM. 

5.3  Electrostatic  Loading  of  Mechanical  Structures 

5.3.1  Doubly  Clamped  Beam 

A  typical  canonical  coupled  electrostatics  and  stress  problem  is  an  elastic  beam  under  an 
electrostatic  load.  The  problem  is  shown  in  Figure  5-7.  The  upper  0.5x10x80  pm  elastic  beam  is 
clamped  at  either  end  and  is  0.7  pm  above  a  rigid  beam.  The  elastic  beam  has  a  Young’s 
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Modulus  of  1.69  xl09Pa,  Poisson’s  Ratio  of  0.3,  and  no  residual  stress.  A  voltage  of  <j>0=  10  V  or 
20  V  is  applied  to  the  upper  beam  and  the  lower  beam  is  grounded  (0  V).  The  electrostatic  model 
calculates  the  electrostatic  pressure  force  on  the  elastic  beam.  The  structural  model  is  used  to 
calculate  the  stresses,  strains,  and  displacements  on  the  elastic  beam.  As  seen  in  Figure  3-1  when 
1 0  V  is  applied  to  the  upper  beam  the  deflection  toward  the  rigid  body  is  small,  with  a  maximum 
displacement  of  0.055  pm  at  beam  center.  With  an  applied  voltage  of  20  V  the  beam  deforms 
onto  the  rigid  body. 


10  nm  80  |xm 

10V  20  V 


Figure  5-7.  Doubly  Clamped  Beam  Under  an  Electrostatic  Load,  with  an  Applied  Voltage  <j)0  = 
10  or  20  V. 

5.3.2  Accelerometer 

Another  problem  shown  in  Figure  5-8(a)  is  an  electrostatically  loaded  plate.  This  is  an  example 
of  an  accelerometer.  A  large  plate  with  an  applied  voltage  of  20  V  is  clamped  by  four  beams. 
The  whole  upper  elastic  structure  has  a  Young’s  Modulus  of  1.69xl09Pa,  Poisson’s  Ratio  of 
0.3,  and  no  residual  stress.  The  upper  plate  is  2.0  pm  above  a  ground  plane.  Figure  5-8(b)  shows 
the  calculated  displacement  contours  on  the  upper  plate  due  to  the  electrostatic  load.  With  20  V 
applied  to  the  upper  plate  a  maximum  deflection  in  the  center  of  the  plate  of  1.83  pm  toward  the 
ground  plane  is  calculated. 


(a)  (b) 

Figure  5-8.  Accelerometer  Under  an  Electrostatic  Load,  (a)  The  geometric  dimensions  and 
problem  set-up.  (b)  The  calculated  displacement  of  the  plate  due  to  the  electrostatic 
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load.  The  displacement  of  the  plate  toward  the  ground  plane  is  maximum  (1.83  pm) 
at  the  center  of  the  upper  plate. 

5.3.3  High  Frequency  Resonator 

Electrostatic  (BEM)/structural  mechanics  coupling  was  demonstrated  on  a  high  frequency 
resonator  shown  below  and  used  for  applications  such  as  high  pass  filters.  A  sinusoidal  driving 
voltage  is  applied  to  a  plate  below  a  resonator  beam.  An  electrostatic  pressure  force  deforms  the 
resonator  beam  as  seen  in  Figure  5-9.  The  deformation  is  coupled  through  a  “coupling  beam”  to 
an  “output  beam”  where  the  change  in  capacitance  between  the  beam  and  the  ground  is  used  to 
detect  its  deformation.  The  contours  show  the  calculated  vertical  displacement  for  one  instance 
in  time. 


Figure  5-9.  Displacement  Field  Contours  on  a  High  Frequency  Resonator 

5.3.4  Electrostatic  Torsional  Micromirror 

Electrostatic  torsional  micro-mirrors  are  employed  to  serve  as  light  modulators  in  devices  such 
as  printers,  scanners  etc.  The  commonly  used  configuration  consists  of  a  flat  plate  made  of 
aluminum  or  polysilicon  that  is  suspended  above  a  ground  plane  by  flexible  hinges.  Here,  an 
aluminum  mirror  of  dimensions  70x40x1.3  pm  (length  x  width  x  thickness)  was  chosen  for 
analysis.  The  torsion  hinges  are  of  dimensions  15x2.5x0.3  pm.  The  plate  is  suspended  2  pm 
above  a  ground  plane  with  embedded  driver  electrodes  that  are  alternately  made  electrically 
active.  A  range  of  voltages  from  0  to  10  V  have  been  analyzed  using  the  FVM/FEM  coupled 
approach.  The  torsion  hinge  was  modeled  as  a  torsion  bar  with  a  specified  spring  constant 
designed  to  describe  the  mechanics  of  the  mirror  as  that  under  a  solid  body  rotation  due  to  an 
applied  moment.  Only  one  half  of  the  mirror  was  simulated  exploiting  the  symmetry  in  the 
problem.  In  the  analysis,  the  standard  far-field  and  symmetry  boundary  conditions  were  specified 
at  appropriate  locations.  The  markedly  disparate  aspect  ratios  in  the  problem  are  an  excellent  test 
of  the  robustness  of  the  employed  approach.  The  coupled  approach  was  able  to  handle  this  issue 
very  well  without  any  stability  or  convergence  problems.  Figure  5-1 0(a)  shows  a  SEM  picture  of 
an  aluminum  micro-mirror  based  upon  which  the  present  geometry  was  roughly  modeled. 
Figure  5- 10(b)  shows  the  mirror,  which  is  grounded,  in  a  deflected  position  due  to  one  of  the  two 
address  electrodes  being  made  active  at  10  V. 
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(a)  micromirror  (SEM  picture) 


(b)  computed  electrostatic  deflection 
Figure  5-10.  Computational  Results  for  an  Electrostatic  Micromirror 


5.3.5  Micromotor  Torque  Curve  Extraction 

The  geometry  and  boundary  element  grid  of  a  micromotor  device  is  shown  in  Figure  5-11.  The 
radius  of  the  rotor  blades  is  50  mm,  and  the  gap  between  rotor  and  stator  is  1 .5  mm.  The  four 
blades  of  the  stator  have  a  potential  of  9.5  V,  the  other  components  are  grounded.  All  parts  are 
conductors.  In  generating  the  grid,  care  was  taken  to  ensure  sufficient  number  of  elements  near 
the  gaps.  The  BEM  solver  was  used  to  compute  the  normal  electric  field  on  the  surface,  which  is 
used  in  turn  to  compute  the  surface  pressure.  The  pressure  is  then  used  to  obtain  the  torque.  The 
computations  were  performed  for  different  rotor  angles.  The  final  torque-angle  relationship  is 
shown  in  Figure  5-12.  The  result  is  in  close  agreement  with  the  Finite  Element  calculations  of 
Barba  et  al.,  1994. 


Figure  5-11.  Geometry  and  Computational  Grid  for  an  Electric  Micromotor 
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Figure  5-12.  Computed  Torque- Angle  Relationship  for  the  Micromotor 

5.3.6  Linear  Lateral  Resonator  Comb  Drive 

The  electrostatic  and  structural  models  were  tested  on  a  linear  lateral  resonator  comb  drive  (Tang 
et  al.,  1992)  shown  in  Figure  5-13.  The  device  has  the  potential  for  many  uses  such  as  an 
accelerometer  or  gyroscope.  A  folded  beam  with  attached  combs  is  placed  between  comb  drives. 
The  comb  drives  have  applied  ac  or  dc  voltages  to  drive  or  sense  the  resonance  of  the  moving 
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folded  beam  structure.  The  sensing  circuitry  takes  advantage  of  the  change  in  capacitance  of  the 
structure  as  it  moves. 

A  calculation  was  performed  on  the  linear  lateral  resonator  comb  drive  shown  in  Figure  5-13. 
For  this  calculation,  one  of  the  comb  drives  is  grounded  and  to  the  other  a  sinusoidal  drive 
voltage  is  applied.  The  ground  plane  and  the  folded  beam  structure  are  also  grounded.  The  folded 
beam  is  only  fixed  at  two  anchor  points  in  the  center  of  the  structure.  The  beam’s  modulus  of 
elasticity  is  1.69  xlO11  Pa  with  a  Poisson’s  Ratio  of  0.3. 

As  the  voltage  increases  at  the  driven  comb  drive  an  electrostatic  pressure  force  pulls  together 
the  comb  drive  and  the  comb  on  the  folded  beam  structure.  The  electrostatic  pressure  force  also 
pulls  the  folded  beam  structure  nearest  it  upwards  away  from  the  ground  plane  due  to  the  electric 
field  asymmetry  introduced  by  the  ground  plane.  This  can  be  seen  in  Figure  5-13  which  shows 
the  structure  at  two  instances  in  time.  The  plotted  contours  represent  the  vertical  displacement 
(displacement  normal  to  the  ground  plane).  The  contours  show  the  resultant  increase  in  tilt  due  to 
the  electric  field  asymmetry  in  the  vertical  direction  as  the  voltage  is  increased. 


comb  drive 
(grounded) 


folded  beam 
(grounded) 


ground 

plane 


comb  drive 
(voltage  applied) 


Figure  5-13.  Linear  Lateral  Resonator  Comb  Drive  with  an  Applied  Sinusoidal  Drive  Voltage  at 
Two  Instances  in  Time 
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Fluid  Damped  Beam  Under  an  Electrostatic  Load 


The  coupled  thermal-fluid-structural-electrostatic  capability  in  CFD-ACE+MEMS  was  tested  on 
a  simple  demonstration  problem.  The  problem  setup  is  shown  in  Figure  5-14.  The  geometry 
consists  of  a  closed  chamber,  5  mm  high  and  70  mm  long,  which  contains  a  compressible  gas  of 
viscosity  equal  to  0.0171  Pa-s  at  an  initial  pressure  P0=  10  Pa  and  temperature  T0=  300  K.  The 
chamber  is  bounded  at  the  top  by  a  flexible  beam  clamped  at  each  edge,  and  by  adiabatic  walls 
on  the  other  three  sides.  The  beam  is  1  mm  thick,  and  has  a  modulus  of  elasticity  of  4000  Pa. 

From  symmetry  considerations,  half  of  the  geometry  is  modeled,  with  symmetry  conditions 
applied  to  the  left  edge. 


sinusoidal  voltage  is  applied  (12000  sin  57tt  V)  to  the  bottom  wall,  and  the  beam  on  top  of  the 
chamber  is  grounded.  The  sinusoidal  varying  electrostatic  force  deflects  the  beam  downward 
compressmg  the  gas.  As  the  gas  is  compressed,  the  pressure  increases,  applying  an  upward  force 
to  the  beam  to  counteract  the  electrostatic  force.  At  each  time  step  the  resulting  deformation  is 
obtained  from  a  balance  of  the  electrostatic,  pressure,  and  inertial  forces. 


The  result  of  the  analysis  was  to  determine  the  transient  response  of  the  fluid  pressure 
emperature,  electrostatic  field,  and  the  deflection  and  stresses  of  the  beam.  To  accomplish  this' 
the  governing  equations  for  the  four  models  (flow,  heat,  electrostatic,  and  structural  analysis)  are 
coupled  and  solved  implicitly  at  each  time  step.  The  time  sequence  of  Figure  5-14  shows  the 
deflection  of  the  beam,  the  gas  flow  as  vectors,  and  the  normal  axial  stress  as  contours  on  the 
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beam.  The  snapshots  were  taken  before,  at,  and  after  the  voltage  maximum.  As  the  beam  moves 
downward  the  gas  is  forced  into  the  clamped  comer  pushing  the  elastic  beam  upward.  The 
energy  equation  was  also  solved  but  the  temperature  increases  were  minimal  (only  0.1  K)  due  to 
the  small  number  of  cycles  simulated. 

A  parametric  study  was  performed  by  varying  the  initial  gas  pressures.  The  graph  of  Figure  5-15 
shows  the  center  beam  displacement  for  three  different  initial  gas  pressures  and  the  lower  wall 
voltage  versus  time.  The  beam  displacement  maximum  is  reduced  with  increasing  initial  gas 
pressure  P0  as  the  increasing  gas  pressure  counteracts  the  electrostatic  pressure  applied  to  the 
beam.  Also,  the  higher  the  initial  pressure,  the  later  the  time  the  beam  reaches  its  maximum 
displacement.  This  is  indicative  of  the  strong  nonlinear  coupling  between  the  different  models. 


P0  =  10  Pa 


Figure  5-14.  Doubly  Clamped  Fluid  Damped  Beam  Under  a  Sinusoidal  Electrostatic  Load. 

Shown  is  a  time  sequence  of  die  beam  displacement,  fluid  velocity  field  (vectors), 
and  normal  axial  stress  (contours  on  beam)  for  the  P0  =  10  Pa  case. 
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Figure  5-15.  Displacement  and  Voltage  Plots  for  a  Doubly  Clamped  Fluid  Damped  Beam  Under 
a  Sinusoidal  Electrostatic  Load.  Shown  are  the  center  beam  displacement  for 
different  initial  gas  pressures  (P0)  and  the  applied  voltage  on  the  lower  chamber 
wall  versus  time. 


5.5  Electrostatic/VOF  Test:  Electrostatic  Extraction  of  Conductive  Fluid  from  Bath 

The  geometry  shown  in  Figure  5.16  was  used  to  demonstrate  the  electrostatic  extraction  of  a 
conductive  fluid  from  a  bath.  The  problem  involves  calculating  gas  flow,  free  surface  flow,  and 
electrostatics.  A  very  large  voltage  is  applied  to  the  upper  wall  of  the  chamber.  The  lower  wall 
and  bath  are  grounded.  The  conductive  fluid  contained  in  the  lower  chamber  contains  charges  on 
its  surface  due  to  the  applied  field.  The  pressure  force  on  the  surface  of  the  conductive  bath  is 
resisted  by  the  surface  tension  of  the  liquid. 
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Figure  5-16.  Time  Evolution  of  the  Electrostatic  Extraction  of  a  Conductive  Fluid  from  a  Bath 

The  free  surface  flows  are  calculated  using  the  volume  of  fluid  (VOF)  method  (Flirt  et  al.,  1981). 
Briefly,  the  VOF  method  solves  an  additional  scalar  quantity  F  which  is  the  volume  fraction  of  a 
component  in  a  cell.  Distribution  of  F  as  a  function  of  time  and  space  is  calculated  and 
determines  the  location  of  the  liquid  cells.  The  main  equation  solved  is 

—  +  u*VF  =  0  (5.10) 

dt 


where  u  is  the  flow  velocity  vector. 

The  VOF  formulation  includes  surface  tension  and  electrostatic  pressure  effects.  To  add  the 
surface  tension  and  electrostatic  pressure  effects,  the  free  surface  of  the  liquid  needs  to  be  located 
and  then  the  forces  applied.  A  conservative  formulation  for  the  surface  tension  forces  calculation 
is  used  (Yang  et  al.,  1998).  It  considers  the  equilibrium  of  a  curved  surface  under  a  surface 
tension  force.  The  net  normal  force  acting  on  the  surface  is  linked  to  the  summation  of  all 
tangential  forces  due  to  surface  tension.  In  the  discrete  form  only  the  boundary  lines  along  the 
surface  are  necessary  to  estimate  the  surface  tension  forces. 

The  results  of  the  calculations  are  shown  in  Figure  5-16  for  different  time  steps.  The  contours  are 
the  electric  potential,  the  vectors  are  the  fluid  flow,  and  the  surface  of  the  liquid  is  represented  by 
the  grey  surface.  As  time  progresses  the  electrostatic  pressure  force  on  the  surface  of  the 
conductive  bath  begins  to  draw  the  liquid  from  the  bath  toward  the  upper  wall  counteracting  the 
surface  tension  forces. 
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5.6  Magnetic  Test  Cases 


5.6.1  Magnetic  Field  Due  to  a  Circular  Current  Carrying  Wire 

An  infinitely  long  wire  of  finite  radius  was  simulated  using  the  magnetic  model.  A  uniform 
current  density  in  the  z-direction  (out  of  the  page)  of  Jz  =  l/p0  A/m2  was  applied  to  the  grey 
region  (a  cross-section  of  an  infinite  wire).  The  calculated  magnetic  field  shown  in  Figure  5-17  is 
rotating  in  the  counter  clockwise  direction  as  one  would  expect  from  the  right-hand-rule.  The 
calculated  and  analytical  values  for  the  z-component  of  the  magnetic  vector  potential  Az  are 
plotted  in  Figure  5-17  (one  plot  is  for  inside  the  wire,  the  other  for  outside  the  wire)  and  show 
very  good  agreement. 


Figure  5-17.  Magnetic  Field  Vectors  Around  Wire  with  Uniform  Current 


Figure  5-18.  Comparison  of  Calculated  (Red)  and  Analytical  (Black)  Solutions  of  Magnetic 
Field  Due  to  an  Infinitely  Long  Current  Carrying  Wire 
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5.6.2  Magnetic  Field  Due  to  Current  Carrying  Bus  Bar  at  a  High/Low  Permeability 
Interface 

As  a  validation  test  case  of  the  magnetic  model’s  multi-material  capability  based  on  the  solution 
of  the  vector  magnetic  potential  equation  a  benchmark  numerical  solution  provided  by  Ida 
(1995)  was  simulated.  This  is  a  2-D  problem  in  which  a  current  carrying  conductor  sets  up  a 
magnetic  field  that  penetrates  into  a  lower  conducting  material,  the  depth  of  penetration  being  a 
function  of  the  relative  permeability  of  the  material.  A  very  fine  grid  of  resolution  300x100  cells 
was  used  to  discretize  the  domain.  The  values  used  in  the  problem  set-up  are  shown  in 
Figure  5-19.  The  contour  plot  of  the  computed  magnetic  vector  potential  is  shown  in 
Figure  5-20.  Comparison  with  the  solution  of  Ida  (1995),  not  shown  here,  shows  excellent 
agreement. 


current  carrying  busbar A/#qm 
'  2  cm  x  1  cm 


free  space 


passive 

conductor 

(n  =  100  (ijj) 


2  cm 


10  cm 


Figure  5-19.  Bus  Bar  Test  Case  Specification 
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Figure  5-20.  Calculated  Vector  Magnetic  Potential  Around  a  Bus  Bar  Between  Materials  of 
Different  Permeability 


5.7  Magnetic  Field  Using  a  Source  Calculated  from  Electric  Model 

Some  testing  was  done  for  the  electric/magnetic  link  when  the  electric  model  solves  the 
conduction  equation  to  supply  the  magnetic  model  with  a  source  current.  This  is  an  important 
link  when  current  sources  of  arbitrary  shape  are  introduced  into  a  problem.  It  would  be 
impossible  for  a  user  to  specify  the  current  source  magnitude  and  direction  in  every  grid  cell  of 
the  source.  The  coupling  of  the  two  models  automates  the  source  current  specification.  Two 
different  source  current  shapes  were  tested:  a  circular  and  square  planar  coil. 


5.7.1  Magnetic  Field  for  Spiral  Conductor 

The  coupling  between  the  electric  and  magnetic  models  is  demonstrated  by  calculating  the 
magnetic  field  around  a  spiral  conductor.  The  terminals  of  the  spiral  are  fixed  at  different 
voltages  and  the  electric  model  calculates  the  current  through  the  conducting  material.  The 
current  is  shown  in  Figure  5-21  as  vectors  colored  with  the  value  of  the  potential.  The  grid 
density  in  the  conductor  is  too  high  to  distinguish  individual  vectors  in  this  plot.  This  current  is 
used  as  a  source  in  the  magnetic  calculation.  The  magnetic  field  vectors  are  shown  in  black  in  a 
plane  perpendicular  to  the  spiral. 
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Figure  5-21.  Magnetic  Field  Due  to  Current  in  Circular  Coil 


The  spiral  has  a  constant  cross  section,  so  the  magnitude  of  the  current  density  is  a  known 
constant  value  everywhere  in  the  conductor.  The  direction  of  the  current,  however,  depends  on 
the  path  length  along  the  spiral. 

5.7.2  Magnetic  Field  Due  to  Planar  Coil 

The  magnetic  actuation  of  structures  is  of  great  interest  due  to  the  larger  distances  which  a 
structure  can  be  moved  when  compared  with  electrostatic  actuation.  The  test  problem,  which  will 
be  used  to  test  this  capability,  is  the  actuation  of  a  plate  due  to  a  magnetic  field  generated  from  a 
planar  coil. 

A  first  step  in  this  problem  is  to  calculate  the  magnetic  field  generated  by  a  planar  coil.  The 
conduction  problem  is  solved  on  the  coil  to  get  a  source  current.  The  source  current  is  then  fed 
into  the  magnetic  solver  which  calculates  the  magnetic  field  due  to  the  current  through  the  planar 
coil.  Shown  in  Figure  5-22  is  the  magnetic  field  due  to  current  flow  in  a  two  turn  planar  coil. 
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Figure  5-22.  Problem  of  Plate  Movement  Due  to  a  Magnetic  Field.  Shown  is  the  (a)  the  electric 
potential  in  the  coil  and  (b)  the  calculated  magnetic  field  generated  from  the  planar 
coil. 


5.8  Magnetic  Actuation  of  a  Beam 

The  magnetic  actuation  of  structures  is  of  great  interest  due  to  the  larger  distances  which  a 
structure  can  be  moved  when  compared  with  electrostatic  actuation.  In  Figure  5-23  is  a 
demonstration  of  the  magnetic  actuation  of  a  beam  in  2-D.  The  deformation  of  the  2-D  beam  was 
due  to  an  induced  magnetic  moment  of  some  attached  magnetizable  material  to  a  beam.  Initially 
a  uniform  magnetic  field  was  applied  to  the  structure.  As  the  material  magnetized,  the  magnetic 
field  near  the  magnetized  region  underwent  variation.  With  a  non-zero  magnetic  moment  and  a 
non-uniform  magnetic  field  in  the  region  of  the  material  a  torque  on  the  material  is  induced.  The 
torque  on  the  material  actuates  the  beam  to  which  the  material  is  attached.  The  displacement  of 
the  beam,  the  magnetic  field  near  the  high  permeability  material,  and  the  torque  on  the  beam  are 
shown  in  Figure  5-23. 


Figure  5-23.  Displacement  Field  Contours,  Magnetic  Field  Lines,  and  Force  Vectors  on  a 
Magnetically  Actuated  Beam 
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5.9  Buoyancy-Driven  Flow  of  a  Conductive  Fluid 


The  Lorentz  force  is  demonstrated  by  the  calculation  of  buoyancy  driven  flow  of  a  conducting 
fluid  in  a  uniform  magnetic  field.  The  bottom  boundary  rotates  with  a  specified  angular  velocity 
and  has  a  uniform  temperature.  The  inner  part  of  the  top  boundary  rotates  in  the  opposite 
direction  of  the  lower  boundary  and  has  a  uniform  temperature  lower  than  that  of  the  bottom 
boundary.  The  outer  part  of  the  top  boundary  is  a  free  surface.  The  temperature  difference 
between  the  boundaries  creates  a  buoyancy-driven  flow  pattern  in  addition  to  the  rotation 
imposed  by  the  walls.  A  uniform  magnetic  field  in  the  opposite  direction  of  gravity  damps  the 
flow  through  the  Lorentz  force.  Figure  5-24  shows  contours  of  the  temperature  and  velocity  in 
the  vertical  direction  with  and  without  the  magnetic  field.  The  velocity  magnitude  is  much 
smaller  when  the  magnetic  field  damps  the  flow. 


Figure  5-24.  (a)  Temperature  and  (b)  Vertical  Velocity  Contours  for  Coupled  Flow/Magnetics 
Solution.  Right  Half  is  without  the  Magnetic  Field,  Left  Half  is  with  the  Magnetic 
Field 


5.10  Electroosmosis 

Electroosmosis  techniques  are  widely  used  for  transport  of  charged  fluids  in  microfluidic 
systems.  Shown  in  Fig.  5-53  is  a  cross  channel  device  used  for  sample  injection  and  separation. 
The  flow  of  material  in  the  device  is  controlled  by  static  fields.  Depending  on  the  voltage  Vs 
applied  to  the  upper  and  lower  reservoir  of  the  cross  channels,  the  flow  of  charged  specie  can  be 
suppressed  (0V),  flowed  into  (37.5V),  or  flowed  out  of  (135V)  the  cross  channels. 
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Figure  5-25.  Electric  Potential  in  Flow  Velocities  for  a  Cross  Channel  Device 
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6.  ANALYSIS  ON  INDUSTRIAL  DEVICES 

The  development  of  the  MEMS  CAD  tool  has  enabled  CFDRC  to  analyze  various  industrial 
problems  and  devices.  The  coupling  of  the  electromagnetic  model  to  the  flow-heat-structure 
model  has  enabled  the  following  simulations. 

6.1  Xerox  Flap  Valve 

Xerox  was  concerned  with  the  voltage  required  to  close  the  flap  on  a  pressurized  valve.  The 
results  are  shown  in  Fig.  6-1  for  four  different  applied  voltages.  The  structural  element  is  a  flap 
over  a  flow  inlet  that  is  fixed  at  the  left,  With  no  flow  or  applied  voltage,  the  flap  lays  over  the 
inlet.  At  OV  with  a  flow,  the  flap  is  deformed,  allowing  flow  through  the  system.  An  electric 
potential  difference  is  applied  between  the  flap  and  the  lower  right  plate.  A  voltage  which  closes 
the  flap,  despite  the  flow,  is  determined. 


300V  325V 


Figure  6-1 .  Flow  Across  a  Flap  with  Different  Applied  Voltages 

6.2  Honeywell  Mesopump 

Honeywell  is  developing  a  mesopump  for  microsystem  applications  (Fig.  6-2).  Specific  pump 
design  parameters  such  as  actuation  voltages,  actuation  times,  and  volumetric  efficiencies  were 
of  primary  interest  to  Honeywell  pump  designers.  The  pump  actuation  is  accomplished  through 
electrostatic  pressure  forces.  This  was  a  flow/structural/electrostatic  problem.  Computational 
results  such  as  valve  closure  times,  membrane  dynamics,  and  pumping  rate  agreed  well  with 
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experiments  conducted  at  Honeywell  (see  Fig.  6-3).  An  example  of  calculated  electric  fields  is 
shown  in  Fig.  6-4. 


Figure  6-2.  Solids  Model  of  a  Mesopump  Cell  with  Boundary  and  Physical  Conditions 


Figure  6-3 .  Comparison  of  Calculated  and  Measured  Actuation  Times  of  Pump  for  Different 
Driving  Voltages 
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Figure  6-4.  Electrostatic  Field  Distribution  in  Pump 

6.3  Honeywell  Cantilever  Beam 

This  test  case  is  an  analysis  of  a  polymer  cantilever  beam.  The  geometry  for  this  device  is 
shown  in  Figure  6-5.  A  kapton  sheet  with  a  thickness  of  25  microns  is  electrostatically  attracted 
toward  a  dielectric  coated  ground  plane.  The  thickness  of  the  dielectric  is  0.2  micron.  In  an 
unstressed  state,  the  kapton  sheet  lies  flat  on  the  dielectric  substrate.  An  upward  directed  force  is 
applied  to  the  kapton  sheet,  over  a  length  denoted  as  L  in  the  figure,  to  raise  it  to  a  specified 
height  h  (typically  20  to  100  microns). 

The  kapton  sheet  is  clamped  on  both  the  right  and  left  edges,  and  thus  nonlinear  structural  effects 
are  important  as  stress-stiffening  effects  become  non-negligible.  In  this  analysis,  the  dielectric 
and  ground  planes  were  assumed  rigid,  and  the  rigid  contact  formulation  was  used  to  model  the 
contact  between  the  kapton  sheet  and  the  dielectric. 

The  purpose  of  the  analysis  is  to  determine  the  force  required  to  raise  the  sheet  for  a  specified 
voltage  and  gap  height.  This  problem  exhibits  an  interesting  nonlinear  effect.  From  a  structural 
standpoint  alone,  more  force  is  required  to  lift  the  sheet  to  a  higher  gap  height,  but  a  smaller  gap 
height  will  require  a  larger  force  to  overcome  the  electrostatic  attraction.  Consequently, 
depending  on  the  structural  properties  and  the  applied  voltages,  there  is  a  possibility  of  a  local 
maxima  in  the  required  force  for  different  gap  heights,  as  opposed  to  a  montonically  increasing 
force  vs.  gap  height  curve. 

Figure  6-6  shows  results  for  gap  heights  of  20  and  60  micron,  with  an  applied  voltage  of  80  V  in 
each  case.  Shown  on  each  plot  are  stress  contours  in  the  kapton  sheet  and  the  electric  field 
through  the  gap,  dielectric,  and  ground  plane.  Figure  6-7  shows  plots  of  the  required  force  per 
unit  area  for  the  case  of  zero  and  80  V.  A  local  maxima  is  seen  in  the  80V  curve  between  a  gap 
height  of  20  and  25  micron.  After  approximately  50  microns  the  electrostatic  attraction  becomes 
negligible  and  the  two  curves  are  nearly  identical. 
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Figure  6-5.  Schematic  of  Electrostatic  Actuator 


Von  Mises  Stress 


Von  Mises  Stress 


Figure  6-6.  Computational  Results  for  h  =  20p  and  h  =  60p 
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Pressure  (N/cmA2) 


Polymer  Electrostatic  Actuator  Analysis 

Required  Pressure  at  0  and  80  V 


Figure  6-7.  Pressure  vs.  Height  Curves  for  0  and  80  Volts 
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7.  IMPLICIT  COUPLING  VIA  MULTI-NEWTON  METHOD 

The  linking  between  physical  models  such  as  the  electric  and  the  structural  models  is  sequential. 
Each  physical  model  is  solved  separately  from  the  others  using  only  boundary  conditions  and 
volume  conditions  the  other  models  calculated  during  their  previous  solution.  When  passing  the 
information  between  models,  linear  relaxation  may  be  applied  to  the  relevant  parameters  to 
improve  convergence.  While  this  linking  is  efficient  for  many  systems,  it  will  fail  when  the 
coupling  becomes  too  stiff.  This  occurs  quite  often  in  solid/fluid  systems  with  small  clearances, 
i.e.,  elastohydrodynamics.  Coupled  electric/structural  systems  with  small  clearances  will  also 
demonstrate  tight  coupling.  In  these  instances,  a  certain  degree  of  implicit  coupling  between  the 
models  will  become  necessary  to  obtain  a  stable,  converged  solution.  One  advanced  implicit 
coupling  methodology  the  “multi-level  Newton  method”  (Aluru  and  White,  1997)  has  been 
implemented  to  ensure  stable  and  robust  convergence.  The  multi-level  Newton  method  has  also 
been  demonstrated  for  electric/structural  coupled  problems. 

7.1  Multi-Newton  Method 

The  implicit  coupling  method  implemented  is  the  “multi-level  Newton  method  being  developed 
by  Jacob  White’s  research  group  at  MIT  (Aluru  and  White,  1997).  The  method  appears  genera 
enough  so  that  it  is  applicable  to  any  physical  model  (structures,  fluids,  heat,  electromagnetics) 
and  any  solution  method  (fvm,  bem,  fern).  The  technique  is  a  black  box  approach  in  that  the 
models  to  be  linked  need  only  have  a  clear  set  of  input  and  output  values  upon  which  t  ey 
mutually  depend.  So,  no  complicated  links  need  be  established  between  the  models. 

In  the  multi-level  Newton  technique,  coupled  equations  are  solved  by  employing  a  Newton- 
Raphson  method.  The  outer  Newton  iteration  solves  a  “Newton-update  linear  system”  which  is  a 
residual  equation  expressing  the  difference  of  the  inputs  and  outputs  of  the  linked  models.  Eac 
Newton  iteration  requires  the  solution  of  the  Newton-update  linear  system,  and  this  linear  system 
can  be  solved  with  a  Krylov-subspace  algorithm  such  as  GMRES.  Each  iteration  of  a  ^J0^' 
subspace  algorithm  requires  a  matrix  vector  product  which  can  be  approximated  using  the  blac 
box  linear  solver  of  each  model.  Since  many  of  these  black  box  solvers  typically  employ 
Newton’s  method,  the  approach  is  referred  to  as  a  multi-level  Newton  method. 

For  instance,  in  an  electro-mechanical  coupling  one  black  box  is  the  structural  mechanical  model 
Rm()  which  inputs  an  electrostatic  force  f  and  outputs  a  grid  displacement  u  =  Rm®-  The  other 
black  box  is  an  electrostatic  model  RE()  which  inputs  a  displaced  grid  u  and  outputs  an 
electrostatic  force  f  =  RE(u).  The  outer  Newton  iteration  solves  a  residual  equation,  that  is  the 

“Newton-update  linear  system”,  of  the  form 


R(u,f) 


Jf-RE(u)| 
\u — R  M  (f )  J 


(7.1) 


The  standard  Newton  technique  is  applied  to  Eq.  (7.1),  in  which  the  residuals  are  expressed  as 
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where  the  ‘  0  ’  represents  the  old  iteration  value. 
The  Jacobian  corresponding  to  Eq.  (7.1)  is  given  by 


J(u,f)  = 


I 

gRM 

8f 


gRE 

du 

I 


which  gives  the  equation 

j£}=-{R(u°-f0)} 


(7.3) 


(7.4) 


to  be  solved  at  each  iteration. 


This  linear  system  can  be  solved  using  only  black-box  solvers  if  Krylov-subspace  based  iterative 
methods  like  GMRES  are  used.  In  these  methods,  the  Jacobian  is  not  explicitly  needed.  Rather, 

V 


matrix-vector  products  such  as  J 
series  is 


are  needed.  The  approximate  form  based  on  the  Taylor 


5f  -  — [RE(u  +  05u )-Rg(u)] 

U 

5U  -  —  [Rm^  +  9§f  )-Riy[(f)] 
0 


(7.5) 


where  0  is  a  perturbation  parameter. 

The  multi-level  Newton  method  may  be  summarized  as  follows: 

1.  Start  with  initial  guesses  for  u,  f  and  calculate  Residual  R(u,  f) 

2.  Solve  Eq.  (7.4)  using  a  Krylov-subspace  iterative  method,  with  the  matrix-vector 
products  approximated  by  Eq.  (7.5) 

3.  Setu  =  u  +  5u;f  =  f+5f 

4.  If  ||5U||  <  su  and  ||5f||  <  Sf,  exit 
Recalculate  new  Residual  R(u,  f)  and  return  to  step  2 


5. 


7.2 


Demonstration  of  Multi-Newton  Method  on  Electric/Structural  Problems 

The  multi-Newton  method  was  used  to  implicitly  couple  the  electric  and  structural  models.  A 
comparison  was  made  between  the  two  coupling  methods:  sequential  relaxation  and  the  multi- 
Newton  method.  The  problem  used  for  comparison  is  a  beam  under  an  electrostatic  load  and 
fixed  at  one  end.  The  beam  is  500  pm  long,  50  pm  wide,  and  14.35  pm  thick  with  a  Young’s 
modulus  of  4.0  xl010N/m2,  a  Possion’s  ratio  of  0.3,  and  no  residual  stress.  A  ground  plate  is 
placed  1  pm  under  the  beam.  The  voltage  applied  to  the  beam  is  varied  to  determine  the  pull-in 
voltage.  The  pull-in  voltage  is  the  voltage  needed  for  the  beam  to  be  pulled  down  and  just  begin 
to  contact  the  ground  plate.  As  voltages  nearer  the  pull-in  voltage  are  applied,  the  coupling 
between  the  electric  and  structural  models  increases. 

Shown  in  Table  7-1  are  results  of  the  comparison.  At  different  applied  voltages,  the  iteration  and 
simulation  time  are  compared  for  the  two  methods.  At  an  applied  voltage  of  7.8  V,  the  Newton 
algorithm  begins  to  out  perform  the  sequential  relaxation  algorithm  (simulation  time  is  less). 


Table  7-1.  Comparison  of  Sequential  Relaxation  and  Multi-Newton  Coupling  Algorithms 


Sequential  Relaxation 

Newton 

V 

DZ 

ITER 

Time 

DZ 

ITER 

Time 

4.0 

-6.0032E-08 

6 

96.85 

-6.0135E-08 

3 

405.94 

5.0 

-9.9109E-08 

8 

128.73 

-9.9507E-08 

3 

422.03 

6.0 

-1.5500E-07 

8 

126.87 

-1.5533E-07 

4 

578.46 

7.0 

-2.4265E-07 

12 

193.16 

-2.4299E-07 

4 

725.35 

7.5 

-3.1792E-07 

20 

319.99 

-3.1907E-07 

5 

795.10 

7.6 

-3.4056E-07 

23 

367.57 

-3.4189E-07 

5 

781.52 

7.8 

-4.2788E-07 

95 

1486.00 

-4.3593E-0 

8 

1291.70 
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8. 


IMAGING,  VISUALIZATION  AND  VALIDATION  OF  MICROSYSTEMS 


Samoff,  in  collaboration  with  the  Exxon  Corporation,  has  established  an  X-ray 
microtomography/microradiography  facility  at  the  Brookhaven  National  Laboratories.  This 
facility  has  the  capability  to  characterize  3-D  microstructure  and  2-D  motion  inside  materials  and 
devices  at  micron  scales.  One  focus  of  this  project  was  to  develop  a  communication  protocol 
between  Samoff  s  imaging  technique  and  CFDRC’s  simulation  tools.  The  developed  protocol 
has  enabled  direct  transfer  of  geometry,  mechanical  motion,  and  velocity  fields  (from 
measurements)  to  be  read  in  directly  to  facilitate  one-to-one  comparisons  between  experiments 
and  simulations. 

8.1  Imaging  and  Visualization 

Imaging  and  visualization  of  microsystems  has  been  carried  out  at  Samoff  Corporation.  The 
initial  investigations  focused  on  the  direct  examination  of  the  3-D  internal  microstructure  and 
microdynamics  in  fluidic  microdevices  using  x-ray  microtomography,  visualization  methods  and 
an  image  processing  supercomputer  (the  Princeton  Engine). 

X-ray  Micro-Imaging:  For  this  project,  Samoff  Corporation  used  the  x-ray  microtomography 
facility  (owned  by  Exxon)  at  the  National  Synchrotron  Light  Source  (NSLS)  at  Brookhaven. 
This  facility  produces  3-dimensional  images  of  the  internal  structure  of  millimeter  sized 
microdevices  with  resolution  approaching  1  micron.  The  measurement  of  the  3-D  microstructure 
and  visualization  of  the  fluid  transport  phenomena  in  microdevices  was  performed. 

Fluid  transport  through  microdevices  can  be  examined  using  digital  subtraction 
microradiography.  This  technique  produces  2-D  projections  of  the  3-D  fluid  flow  through  the 
channels  where  the  background  structure  of  the  specimen  has  been  digitally  removed  from  the 
image.  Liquids  such  as  water  are  injected  into  the  channels  by  a  syringe  pump  operating  at  a 
constant  flow  rate  ranging  from  pl/min  to  ml/min.  A  contrasting  agent  such  as  chlorine  or  iodine 
is  added  to  one  of  the  liquids  since  at  20keV  there  is  little  contrast  between  the  hydrocarbons  and 
water  in  micron  sized  paths.  The  experiment  proceeds  by  taking  “pictures”  (radiographs)  every 
few  seconds  (eventually  at  high  speed)  until  enough  frames  are  collected  to  visualize  the 
displacement.  Typically,  over  1000  frames  are  accumulated.  Several  frames  of  the  dry  specimen 
are  collected,  integrated  and  then  subtracted  from  the  radiographs.  The  results  are  displayed  on  a 
monitor  and  eventually  recorded  on  videotape. 

These  flow  experiments  are  coupled  with  a  measurement  of  the  static  3-D  microstructure  of  the 
microdevice  using  x-ray  microtomography.  This  technique  provides  complementary,  non¬ 
destructive,  morphological  characterization  of  the  device  and  the  fluids  inside. 

Supercomputing  and  Visualization:  The  large  quantity  and  size  of  the  data  sets  generated 
dictate  the  use  of  an  image  processing  supercomputer  for  the  automatic  analysis  and  display  of 
information.  In  particular,  Samoff  Corporation  has  developed  the  Princeton  Engine  (PE),  a 
massively  parallel  computer.  It  is  designed  specifically  for  real-time  image  processing  and 
applications  requiring  extremely  high  input  bandwidth.  The  PE  is  used  for  real-time,  interactive 
volume  rendering  of  the  microtomographic  volume  data  sets.  Flow  characteristics  can  be  better 
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understood  when  viewed  in  the  context  of  the  entire  three-dimensional  microstructure. 
Techniques  for  visualizing  the  structure  volumetrically  include  real-time  interactive  rotation.  For 
example,  a  1283  volume  data  set  can  be  interactively  rotated  at  a  rate  of  30  frames/sec.  It  is  only 
when  these  data  sets  are  set  in  motion  at  sufficiently  high  frame  rates,  that  the  3-D  geometry  of 
the  microstructure  becomes  apparent.  In  addition,  features  of  interest  can  be  exposed  within  a 
volumetric  data  set  by  manipulating  the  opacity  and  color  maps. 

8.2  Three-dimensional  Geometry  Extraction 

The  3-D  microstructure  data  and  the  2-D  motion  data  are  obtained  using  two  types  of 
experiments  described  above,  both  of  them  non-destructive.  Microtomography  experiments 
involve  the  determination  of  the  actual  structure  of  a  fabricated  MEMS  device.  Since  the  field  of 
MEMS  fabrication  is  relatively  new,  and  the  dimensions  are  so  small,  it  is  unclear  how  closely  a 
manufactured  device  matches  the  original  design  specification.  In  this  process,  images  of  a 
sample  are  taken  at  fixed  intervals,  as  the  stage  upon  which  the  sample  sits  is  rotated.  These 
images  are  then  used  as  input  to  a  numerical  technique  called  reconstruction,  which  produces  a 
single  3-D  image  of  the  sample  with  resolution  approaching  1  micron.  This  image  can  then  be 
used  as  the  input  for  a  computational  analysis  program  to  ensure  that  the  actual  geometry,  rather 
than  an  idealization,  is  used. 

Three  different  methods  of  geometric  extraction  were  compared.  Each  method  is  described 
briefly  below. 

Method  1.  Isosurfacing  via  Marching  Cubes 

Isosurfacing  involves  the  construction  of  a  connected  set  of  polygons  at  a  specified  threshold 
value  of  the  image  function.  Essentially,  we  are  choosing  to  connect  neighboring  voxels  whose 
value  equals  the  threshold  value.  The  threshold  value  is  chosen  such  that  it  isolates  a  single 
component  in  the  data  set.  The  basic  algorithm  steps  from  one  end  of  the  data  set  to  the  other, 
“marching”  through  each  voxel.  Only  information  local  to  the  current  voxel  is  used  in  the  surface 
generation.  As  a  result  of  this,  isosurfacing  tends  to  extract  imaging  noise  and  can  suffer  from  the 
formation  of  false  holes  and  spurious  artifacts.  Figure  8-1  illustrates  this  technique  applied  to  a 
channel  comprised  of  an  asymmetric  expansion.  Note  the  apparent  roughness  of  the  surface  and 
the  presence  of  small  outlier  surfaces  surrounding  the  main  channel.  The  surface  roughness  is  not 
’’real”,  but  rather  stems  from  the  fact  that  the  marching  cubes  algorithm  fits  the  surface  to  the 
noise  in  the  underlying  data 
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Figure  8-1.  Extracted  Geometry  Algorithm:  Isosurfacing  via  Marching  Cubes 

An  order  of  magnitude  increase  in  speed  can  be  achieved  by  applying  the  Marching  Cubes 
algorithm  to  only  the  segmented  regions  and  ignoring  all  other  points.  The  accuracy  of  the 
algorithm  was  enhanced  during  the  project  by  modifying  the  Geometry  Extraction  algorithm  to 
utilize  correctly  the  original  tomographic  gray-scale  information  (rather  than  arbitrary  gray  scale 
values  as  was  done  previously)  in  the  vicinity  of  the  segmented  object  surface.  This  modification 
has  resulted  in  reduced  artifacts  and  smoother  triangulated  surfaces.  Figure  8.2  compares  the 
results  from  the  old  and  new  algorithms  for  the  Redwood  Microsystems  Normally  Open 
Microvalve. 


(a) 


(b) 


Figure  8-2.  Comparison  of  Old  and  New  Marching  Cubes  Algorithm  on  the  Redwood  Systems 
Normally  Open  Microvalve 
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Method  2.  Active  Contours 

In  contrast  to  isosurfacing,  active  contour  algorithms  seek  to  fit  piecewise  low-order  polynomials 
to  feature  boundaries  (edges)  in  the  image  function.  Two  initial  steps  are  performed:  the  user 
specifies  initial  estimates  of  the  feature  boundaries  and  the  input  image  function  is  filtered  with 
an  edge-enhancing  operation.  Then,  via  an  iterative  process,  the  initial  polynomial  estimates  are 
shifted  such  that  the  average  proximity  of  the  polynomial  to  the  edge  features  is  minimized.  If 
the  initial  estimates  are  not  separated  from  the  desired  features  by  any  local  artifacts,  the  fitted 
polynomials  will  appear  to  “cling”  to  feature  boundaries.  This  process  can  be  performed  on  each 
slice  of  the  volume  and  the  results  meshed  together.  Figure  8-3  shows  this  technique  applied  to 
the  channel  data  set.  We  have  found  the  active  contour  approach  to  work  very  well  for  a  certain 
class  of  imagery,  generally  characterized  by  sharply  defined  edges  that  form  a  closed  surface. 
However,  in  cases  where  the  geometry  is  complex  or  the  definition  between  the  desired  feature 
and  its  background  is  not  clear,  this  method  breaks  down.  Additionally,  there  is  not  a  great  deal 
of  sensitivity  to  the  initial  polynomial  estimate,  making  this  technique  not  easily  amenable  to 

automation. 


Figure  8-3.  Extracted  Geometry  Using  an  Active  Contours  Algorithm 

Method  3.  Object  Segmentation 

Object  segmentation  involves  the  application  of  a  series  of  filters  to  the  input  volume  to  produce 
an  intermediate  “binary”  volume.  This  intermediate  volume  represents  basic  structures  in  the 
input  in  a  form  in  which  the  effects  of  imaging  noise  have  been  reduced.  The  intermediate 
volume  can  then  be  fitted  with  geometrical  primitives  via  isosurfacing.  The  initial  output  is  well 
behaved  in  terms  of  its  sensitivity  to  noise  and  requires  no  initial  feature  estimates.  The  filtering 
operations  include  an  initial  thresholding  operation  to  convert  the  input  image  into  a  binary 
image,  followed  by  a  connected  component  analysis.  Here,  the  idea  is  to  isolate  significant 
features  and  avoid  the  outlier  problems  noted  above  with  pure  isosurfacing.  After  determining 
the  features  of  interest,  several  morphological  operations  are  applied  which  have  the  effect  of 
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filling  holes  created  by  digitization  noise  and  removing  thin  isthmus’  which  falsely  connect 
larger  features.  Finally,  a  smoothing  operation  is  performed  which  ameliorates  the  effects  of 
imaging  noise  on  the  overall  surface.  Figure  8-4  illustrates  object  segmentation  applied  to  the 
channel  data  set.  In  terms  of  general  usability,  we  have  found  the  object  segmentation  technique 
to  be  superior  to  simple  isosurfacing  and  active  contours.  However,  no  one  single  technique 
works  well  in  all  cases  and  it  is  useful  to  have  a  number  of  algorithms  from  which  to  choose 
based  on  the  class  of  imagery  to  be  analyzed. 


Figure  8-4.  Extracted  Geometry  Using  an  Object  Segmentation  Algorithm 

Additional  Constraints 

As  explained  above,  extracting  the  flow  channel  surface  from  tomography  data  is  basically  a 
process  of  finding  a  set  of  polygons  which  best  fit  prominent  features  in  the  data.  However,  since 
this  data  will  be  used  as  the  geometry  input  to  a  CFD  solver,  additional  constraints  must  be 
added  to  this  process.  One  constraint  is  that  any  degenerate  triangles  must  be  removed  from  the 
data.  The  other  constraint  is  that  the  surface  set  should  reflect  the  natural  boundary  conditions  of 
the  system,  e.g.  planar  inlet  and  outlet  regions.  Quite  a  bit  of  the  work  during  this  project 
involved  imposing  these  additional  constraints  on  the  tomography  data  so  that  they  could  be  used 
by  CFDRC’s  flow  solvers.  Samoff  satisfied  these  constraints  by  employing  a  mesh  decimator 
and  other  filters  to  the  original  triangulated  surface,  until  the  degenerate  triangles  were  removed 
and  the  inlet  and  outlet  surfaces  were  planar. 

Geometric  Accuracy 

After  the  geometries  were  extracted,  a  method  was  developed  to  give  a  quantitative  estimate  of 
the  accuracy  of  the  process.  The  approach  here  was  to  create  a  “ground-truth”  data  set,  for  which 
the  exact  locations  of  spatial  features  are  known.  Since  the  physical  process  (i.e.  transmission 
radiography)  employed  in  the  generation  of  the  imagery  is  known,  the  effect  can  be  modeled  in 
software  and  an  “ideal”  radiograph  of  some  geometry  can  be  calculated.  This  synthetic  data  set 
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can  then  be  employed  as  the  input  to  the  reconstruction  algorithm,  which  converts  the  sequence 
of  images  into  a  3-D  volume  data  set.  The  geometry  extraction  algorithms  are  then  applied  to  the 
synthetic  volume  and  the  resulting  surface  descriptions  are  compared  with  those  from  the 
initially  specified  geometry.  From  this  comparison,  conclusions  can  be  drawn  regarding  the 
effectiveness  of  the  reconstruction  and  geometry  extraction  processes.  Figure  8-5  shows  a 
channel  data  set  of  an  original  model  and  the  resulting  tomography. 


Figure  8-5.  Comparison  of  the  Geometry  Extracted  from  the  Synthetic  Tomography  (Blue) 
versus  the  Original  CAD  Model  (Yellow).  The  agreement  is  very  good  with  the 
exception  of  a  slight  rounding  of  the  edges. 

From  analysis  of  these  results,  the  following  points  are  evident: 

The  comers  of  the  channel  walls,  which  are  right  angles  in  the  model,  are  slightly  rounded.  This 
rounding  is  about  5%  of  the  length  scale,  and  is  a  result  of  a  smoothing  operation  performed 
during  geometry  extraction. 

The  flat  areas  of  the  channel  walls  are  represented  as  perfectly  flat  surfaces  in  the  extracted 
geometry.  This  is  in  contrast  to  slight  but  consistent  roughness  seen  in  real  data.  We  conclude 
that  geometry  extraction  does  not  appear  to  contribute  significantly  to  error  in  spatially  uniform 
areas. 

The  areas  of  gentle  curvature  are  recovered  to  first  order,  but  not  exactly.  The  deviation  is  of  the 
same  order  as  that  noted  for  the  comer  areas,  but  appears  to  be  more  random  in  nature.  This 
effect  may  be  a  result  of  the  representation  of  smooth  curves  by  discrete  polygons  in  the  CAD 
model.  However,  this  uncertainty  is  still  less  than  that  attributable  to  imaging  noise. 

In  summary,  it  appears  that  the  uncertainty  generated  by  the  reconstruction  and  geometry 
extraction  processes  are  small  in  relation  to  the  experimental  uncertainty  of  the  image  acquisition 
hardware. 
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8.3  Two-dimensional  Flow  Field  Extraction 


In  microradiography  experiments,  the  object  in  question  is  stationary  and  X-ray  images  are  taken 
at  fixed  temporal  intervals  until  enough  frames  are  collected  to  visualize  some  dynamic  process 
(e.g.  actuation  of  a  flap  valve,  fluid  motion,  etc.).  In  studying  fluid  flow,  a  neutrally  buoyant 
emulsion  is  used  to  track  fluid  motion.  The  recovery  of  the  flow  field  from  these  images  is  a  non¬ 
trivial  matter  because  of  local  ambiguities  in  the  apparent  displacement  of  image  intensities  and 
the  effects  of  noise.  Samoff  constrains  the  flow  recovery  process  by  making  use  of  physical  and 
geometric  constraints  derived  from  fluid  mechanics. 

Currently,  three  constraints  are  imposed  on  the  flow  recovery  process.  The  first  constraint 
enforces  local  rigidity  of  the  recovered  flow  (i.e.  infinitesimal  material  particles  only  translate 
and  rotate).  The  second  constraint  enforces  local  smoothness  of  the  recovered  flow  (effectively 
constraining  the  velocity  change  from  point  to  point).  This  constraint  is  particularly  important  in 
minimizing  the  effects  of  image  noise.  The  third  constraint  is  related  to  boundary  conditions  on 
the  recovered  flow.  Specifically,  the  flow  must  stay  in  its  channel  and  the  velocity  must  go  to 
zero  at  the  channel  surface. 

Variational  calculus  is  used,  along  with  the  three  constraints  mentioned  above,  to  yield  partial 
differential  equations  which  relate  image  density  to  the  flow  field.  Solution  of  the  resulting 
equations  results  in  a  2-D  flow  field,  representing  the  averaged  flow  across  the  height  of  the 
device.  These  equations  are  solved  numerically  on  a  grid  that  can  be  composed  of  arbitrary 
collections  of  polygons  to  approximate  complex  device  geometries. 

8.4  Comparison  of  Measured  and  Calculated  Flow  Data 

The  image  data  from  Samoff  is  provided  as  a  triangulated  surface  in  STL  (StereoLithography) 
format.  This  data  then  provides  the  geometric  input  to  the  flow  solver.  The  triangulated  surface 
can  either  be  used  as  the  surface  grid  for  the  CFD  mesh,  or  may  be  used  to  define  the  surface 
geometry  onto  which  a  separate  grid  is  imposed.  In  the  early  stages  of  this  project,  the  CFDRC 
geometry  and  grid  generation  software  was  not  able  to  use  STL  files  as  a  surface  geometry 
definition,  and  thus  we  attempted  to  use  the  STL  files  as  a  triangulated  surface  mesh  from  which 
a  volume  mesh  was  derived.  This  process  was  successful,  but  the  results  obtained  from  the  flow 
solver  were  not  satisfactory.  The  surface  grid  obtained  from  Samoff  used  triangles  which  were 
all  nearly  the  same  size.  This  is  a  good  approach  for  geometry  specification,  but  it  is  a  poor 
method  for  constructing  a  surface  mesh,  since  there  is  not  a  finer  grid  in  regions  of  the  domain 
where  large  flow  parameter  gradients  may  occur.  The  result  of  this  method  was  a  grid  which  had 
a  fairly  large  number  of  grid  cells,  but  was  not  fine  enough  in  certain  regions  to  capture 
important  flow  details.  Basically,  we  were  asking  more  of  the  tomography  reduction  procedure 
than  it  was  intended  for.  As  a  result,  for  the  remainder  of  the  cases  we  constructed  the  grid  from 
the  sketches  provided  by  Samoff. 

Figure  8-6  shows  four  flow  channels  for  which  image  data  was  provided  to  CFDRC  from 
Samoff.  These  flow  channels  each  were  200  micron  high,  and  were  etched  in  a  silicon  substrate 
and  covered  with  glass.  The  tomography  images  for  these  geometries  are  shown  in  Figure  8-7. 
The  computational  grids  used  in  for  the  CFD  runs  are  given  in  Figure  8-8.  These  grids  were 
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created  “from  scratch”  starting  with  the  sketches  shown  in 
shown  in  Figure  8-7. 


(a) 


(c) 


(d) 


Figure  8-7.  Visualization  of  the  STL  Files  of  the  Flow  Channels  Received  from  Samoff 
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Channel  12  CFD  Grid 


(C) 


(d) 


Figure  8-8.  Computational  Grids  for  the  Channel  Geometries 


After  the  CFD  solution,  the  3-D  flow  fields  must  be  averaged  across  the  channel  height  for 
comparison  with  the  Samoff  experimental  data.  A  special  post-processing  program  was  written 
to  perform  this  task.  In  addition,  another  program  was  developed  to  convert  the  flow  data 
provided  by  Samoff  into  a  format  which  could  be  read  by  CFDRC’s  post-processing  graphics 
program,  CFD-VIEW.  Figure  8-9  through  Figure  8-11  shows  the  comparison  of  the  averaged 
flow  fields  from  the  Samoff  experiments  and  the  CFDRC  simulations  for  three  of  the  channels. 
Flow  data  was  not  available  for  channel  26  so  Figure  8-12  shows  only  the  CFD  predicted 
velocity  and  stresses.  For  each  case,  velocity  vectors  are  shown  colored  with  the  magnitude  of 
the  predominant  velocity  component.  As  seen,  the  comparison  is  very  good,  both  qualitatively 
and  quantitatively. 
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Channel  12  Experimental  Results 


Channel  12  Numerical  Results 
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Figure  8-9.  Numerical  and  Experimental  Velocity  Fields  for  Channel  12 


Channel  1 4  Experimental  Data 


Channel  14  Numerical  Data 
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Figure  8-10.  Numerical  and  Experimental  Velocity  Fields  for  Channel  14 


Channel  21  Experimental  Data 


Channel  21  Numerical  Data 


Figure  8-11.  Numerical  and  Experimental  Velocity  Fields  for  Channel  21 
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Channel  26  Numerical  Data 


Figure  8-12.  Predicted  Velocity  Field  for  Channel  26 

During  the  course  of  this  project,  work  was  progressing  on  a  different  project  within  CFDRC  for 
the  development  of  an  automatic  grid  generation  program  for  external  flow  aerodynamics.  After 
the  initial  development  work,  the  program  capabilities  were  extended  to  accommodate  internal 
flows.  The  program  developed  from  this  work,  CFD-CARTESIAN,  can  create  adaptive 
Cartesian  meshes  from  surface  triangulations  in  various  formats,  including  STL  format,  ft  uses 
an  omni-tree  data  structure  to  support  anisotropic  grid  adaptation,  and  also  has  the  ability  to 
create  viscous  sublayers  for  resolution  of  boundary  layers. 

To  demonstrate  the  use  of  this  program,  Channels  12  and  14  were  gridded  with  CFD- 
CARTESIAN.  These  grids  are  shown  in  Figure  8-13.  For  channel  14,  the  mesh  was  used  in  a 
flow  simulation  using  CFD-ACE+.  Figure  8-14  shows  the  flow  results  at  the  center  plane  (these 
results  were  not  averaged  across  the  thickness  because  the  post  processing  program  which  does 
this  does  not  accommodate  adaptive  Cartesian  meshes).  The  predicted  maximum  velocity  is  1.4 
times  the  average  velocity  shown  in  Figure  8-10  (compared  to  the  theoretical  value  of  1.5). 
These  results  compare  favorably  considering  the  difference  in  grid  density  between  the  two  runs. 


Figure  8-13.  CFD-CARTESIAN  Calculated  Conversion  of  STL  Format  Files  into  Computational 
Grids  for  Channel  12  and  14 
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Channel  14  -  Adaptive  Cartesian  Mesh 
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Figure  8-14.  Resulting  CFD  Calculation  on  Channel  14  Using  the  Computational  Grid  Generated 
by  CFD-C ARTESIAN  from  the  STL  Formatted  File 


8.5  Comparison  of  Measured  and  Calculated  Electrostatic  Deformation 

A  search  was  made  to  find  a  sufficiently  experimentally  characterized  electrostatically  actuated 
device  for  which  the  developed  CFD-ACE+MEMS  software  could  be  validated.  Ideally,  the 
experimental  characterization  would  contain:  (i)  measurement  of  the  fabricated  device  geometry; 
(ii)  a  specification  of  the  idealized  geometry  of  the  device  such  as  specified  by  a  fabrication 
definition  file  such  as  the  “CIF”  format;  (iii)  the  electronic  drive  and  sensing  circuitry  the  device 
was  operated  in  during  testing  and  visualization;  and  (iv)  measurements  of  the  motion  of  the 
device  such  as  obtained  by  analysis  of  experimental  imagery. 

The  data  Samoff  can  supply  includes:  (i)  images  of  devices  in  operation  at  eight  intervals  along 
the  sinusoidal  drive  voltage  for  different  drive  frequencies;  (ii)  (x,y,t)  motion  data  for  the  eight 
intervals  as  extracted  from  the  2-D  images;  (iii)  CIF  format  files  of  the  manufacturing 
specification  of  the  device;  (iv)  input/output  signal  plot  over  a  range  of  frequencies  to  determine 
Q;  and  (v)  the  z  component  of  motion  data  (in  addition  to  x,y)  at  resonance  although  guessed  at 
by  analyzing  2-D  images. 

Two  such  data  sets  of  previously  visualized  devices  were  considered  for  this  validation  task.  One 
was  a  linear  lateral  resonator  comb  drive  with  the  imaged  field  of  view  of  the  device  shown  in 
Figure  8-1 5(a)  and  the  other  was  an  angular  comb  drive  (manufactured  by  Exponent  Failure 
Analysis  Corporation)  whose  imaged  field  of  view  is  shown  in  Figure  8-1 5(b). 
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Figure  8-15.  Representative  Visualization  Images  of  Two  Electrostatically  Actuated  Comb 
Drives,  (a)  Linear  Lateral  Resonator  Comb  Drive  and  (b)  Angular  Comb  Drive 

There  was  not  enough  imaging  information  on  the  linear  lateral  resonator  comb  drive  (Figure 
8- 15(a))  to  simulate  the  device.  The  imaged  region  of  the  device  has  clipped  off  the  folded 
beams  on  either  side  of  the  comb  structures.  The  length  of  these  folded  beams  is  an  essential 
parameter  which  controls  the  resonant  states  of  the  device.  The  “CIF”  file,  which  establishes  the 
processing  specification  used  to  manufacture  the  device,  was  unavailable  so  even  an  idealized 
geometry  of  this  device  was  unattainable. 

The  other  device  under  consideration  is  an  angular  comb  drive  (Figure  8-15(b)).  Again  the 
imaged  data  clips  off  some  essential  geometric  features  of  the  device.  One  is  the  length  of  the 
beam  in  the  lower  right-hand  side  of  the  image  which  is  a  deforming  beam  whose  length  will 
have  a  large  impact  on  device  operation.  Another  piece  of  geometry  missing  is  the  length  of  the 
fingers  of  the  comb  drive  toward  the  rotational  center  of  the  device.  The  “CIF”  file  of  this  device 
is  available  however  (Figure  8-16)  so  that  an  idealized  geometry  of  the  device  can  be  generated 
for  the  numerical  simulation. 


84 


Figure  8-16.  Visualization  of  CIF  File  Specifying  MUMPS  (Multi-User  MEMS  Processes) 
Manufacturing  Process  for  Construction  of  Angular  Resonator 


Figure  8-17.  Angular  Resonator  Grid  Used  for  Simulation 
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The  grid  for  the  angular  resonator  (Figure  8-17)  was  generated  by  reading  in  the  CIF  formatted 
file  into  CFD-GEOM.  GEOM  converted  the  file  into  lines  and  points.  All  grid  specifications  had 
to  be  done  by  hand.  A  structured  grid  was  used  near  the  interlaced  comb  fingers  for  which  BEM 
is  more  accurate.  An  unstructured  grid  was  used  for  the  center  resonating  arm.  Prisms  were  used 
for  the  unstructured  volume  mesh. 


A  sinusoidal  voltage  Vmsin(27if  t  +  <j>)  with  Vm=35.15V,  f=  20387.9 Hz,  and  <j>  =  0  was  applied 
to  the  fixed  comb  on  the  right  in  Figure  8-17.  All  other  structures  were  grounded.  In  the 
structural  model  the  faces  of  the  center  resonating  arm  connected  to  the  contact  pad  were  given  a 
fixed  displacement.  The  structural  mechanics  model  was  solved  only  on  the  center  resonating 
arm  volume  grid  (13315  prisms).  The  electric  BEM  model  was  solved  only  for  the  structured 
portion  of  the  ground  plane  under  the  right  fixed  comb  and  the  fingers  and  connecting  faces  of 


t  =  T/6 


t  =  T/4 


Figure  8-18.  Calculated  Deformation  of  the  Angular  Resonator  at  Two  Different  Times.  Color 
Contours  are  Displacement  Magnitudes. 
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Shown  in  Figure  8-18  is  the  calculated  deformation  of  the  angular  resonator  at  two  different  . 
times.  The  color  contours  are  displacement  magnitudes  with  a  peak  of  5.5  pm  occurring  at  time 
T/4  where  T  is  the  period  (1/f). 

The  calculated  data  cannot  be  reliably  compared  to  the  measured  data  from  Samoff.  The  actual 
Young’s  Modulus  of  the  fabricated  device  is  not  known.  The  actual  geometry  of  the 
manufactured  device  was  not  gridded,  only  the  idealized  geometry  from  an  electronic  file  was 
used.  Also,  if  any  detecting  voltage  was  applied  to  the  left  comb,  it  would  effect  the 
displacement  of  the  drive.  The  presence  or  absence  of  the  detecting  voltage  on  the  left  fixed 
comb  is  not  known. 

The  experimentally  measured  displacement  is  calculated  by  interpolating  changes  in  the  imaged 
device.  The  experimental  determined  displacement  of  the  point  in  the  center  of  the  red  box  in 
Figure  8-15(b)  is  in  peak-to-peak  values:  Ax  =  666  nm  peak-to-peak.  Ay  =  2090  nm  peak-to- 
peak,  and  Az  =  55  nm  peak-to-peak.  The  calculated  values  using  the  CAD  tool  are:  Ax  =  1621 
nm,  Ay  =  4955  nm,  and  Az  =  1  nm.  The  agreement  between  these  numbers  is  not  very  good  but 
could  be  attributes  to  many  things  as  noted  above. 
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9.  CONCLUSIONS  AND  PLANS  FOR  FUTURE  WORK 

This  chapter  presents  the  conclusions  from  the  current  project  and  the  plans  for  future  extensions 
and  applications  of  the  models. 

9.1  Summary  of  Accomplishments 

The  completion  of  this  project  has  resulted  in  the  first-ever  integration  of  fluidic,  thermal, 
chemical,  electromagnetic  and  mechanical  models  within  a  single  software  environment,  CFD- 
ACE+MEMS.  In  addition  to  the  physical  models,  CFD-ACE+MEMS  has  geometry/grid 
generation  modules,  scientific  data  visualization  module  and  a  fluid/matenal  properties  database. 
Depending  on  the  problem  being  solved,  the  user  has  the  option  to  either  have  a  loose  coupling 
between  the  physical  models  or  a  tight,  implicit  coupling  (necessary  for  highly  stiff  systems)^ 
All  of  the  software  modules  have  been  developed  by  CFDRC  and  will  be  marketed  an 

supported  by  CFDRC. 

In  addition  to  development,  the  models  and  the  coupling  between  the  modLiles  were  tested  and 
verified  on  several  simple  as  well  as  complex  device  systems.  The  FVM  and  BEM  models  for 
electrostatics  were  benchmarked  for  different  applications.  This  project  has  also  produced  a 
communication  protocol  between  Samoff  s  unique  microtomography  technology  (for  visualizing 
device  microstructure,  motion  and  velocity  fields)  and  CFDRC  s  simulation  so  tware.  is 
protocol  enables  the  generation  of  geometry  and  computational  grid  directly  from  the  imaged 
data  and  facilitates  one-to-one  comparison  between  measured  and  computed  device  behavior. 
The  agreement  between  model  and  measurements  was  evaluated  for  several  microfluidic  devices 

as  well  as  an  electrostatic  rotary  comb  drive. 

The  papers  published  resulting  from  this  effort  are  Stout  et  ah,  1999;  Athavale  et  al,  1999;  and 
Krishnan  A.,  1999. 

Finally,  CFDRC  established  working  relationships  with  several  MEMS  device  manufacturers 
(during  this  project)  to  demonstrate  the  simulation  tools  for  design  of  industrial  devices.  ese 
collaborations  are  summarized  below. 

9.2  Industrial  Collaborations  and  Software  License 

This  project  has  enabled  CFDRC  to  perform  several  recent  modeling  studies  for  MEMS 
companies  using  the  coupled  electrostatic  capability.  These  are  summarized  below: 

Micropump  for  Honeywell:  CFDRC  performed  3-D  simulations  of  coupled  fluid- 
structural-electrostatic  phenomena  in  a  micropump  currently  being  prototyped  by 
Honeywell.  These  simulations  were  successful  in  highlighting  some  of  the  complexities 
of  the  interaction  between  the  above  phenomena  and  enabled  Honeywell  designers  to 
improve  the  performance  of  the  micropump. 

Electrostatically  Activated  Beam  for  Honeywell:  Simulations  were  performed  for 
Honeywell  to  analyze  the  deflection  of  a  cantilever  beam  for  different  actuation  voltages. 


88 


Non-linear  effects  were  also  simulated.  Mechanical  forces  (as  a  function  of  gap  height) 
for  different  voltages  were  computed. 

•  Microvalve  for  Xerox:  CFDRC  performed  simulations  of  coupled  fluid-structural- 
electrostatic  interactions  in  Xerox’s  microvalve.  The  simulations  clearly  showed  the 
highly  non-linear  nature  of  the  valve  operation  and  also  gave  insights  as  to  the  extent  of 
voltage  drop  necessary  for  optimal  performance  of  the  valve. 

•  Electrophoresis  Applications:  CFDRC  is  currently  working  with  Oak  Ridge  National 
Lab.,  Caliper  Technologies,  Aclara  Biosciences,  etc.  to  demonstrate  CFDRC’s 
electrophoresis  model  that  accounts  for  interactions  between  charged  species  transport 
and  the  applied  electrostatic  field. 

Several  organizations  such  as  Redwood  Microsystems,  Lucas  Novasensor,  YSI,  Motorola,  Oak 
Ridge  National  Laboratory,  BioDOT,  Stanford  University,  University  of  California  (Berkeley), 
University  of  Washington,  etc.  are  currently  using  CFD-ACE+MEMS  for  microsystem 
applications.  Additionally,  several  microelectronics  companies  such  as  Wacker-Siltronic, 
Aixtron,  Motorola,  Applied  Materials  and  Novellus  have  purchased  software  licenses  from 
CFDRC  for  coupled  fluid-thermal-structural-electromagnetic  simulations.  Their  applications 
include  electromagnetic  induction  heating,  Lorentz  force  stabilization  of  melt  flow  in  crystal 
growth,  low-pressure  plasma  transport  and  electroplating.  CFDRC  expects  a  large  market  for  this 
coupled  simulation  capability.  CFDRC’s  software  is  one-of-a-kind  in  its  ability  to  closely 
integrate  these  models. 

9.3  Future  Work 


This  project  has  laid  the  foundation  for  CFDRC  to  extend  its  capabilities  to  (i)  the  development 
of  reduced/parametric  models,  and  (ii)  simulation  of  systems  of  microdevices  using  the  reduced 
models.  Two  key  projects  (currently  ongoing  at  CFDRC)  that  will  build  upon  this  foundation  are 
as  follows  : 

•  DARPA  BAA  97-17:  “Generation  of  Reduced  Parametric  Models  of  Microdevices  from 
High  Fidelity  Tools  for  System  Level  Composite  CAD”  focuses  on  building  a  high-fidelity 
simulation  environment  for  generating  reduced  models  for  MEMS  devices  and  systems.  This 
tool  will  perform  mixed-dimensionality  simulations  (i.e.,  simultaneously  simulate  a  system 
of  a  micro-channel  and  a  micro-mixer  using  a  0-D  model  for  the  channel  and  a  2-D  model  for 
the  mixer).  Mixed-dimensionality  simulations  will  be  performed  to  determine  dynamic 
responses  of  devices  which  will  then  be  used  to  derive  reduced  order  ‘lumped-parameter’ 
models. 

•  DARPA  BAA  97-39:  “Mixed-Dimensionality  VLSI-Type  Configurable  Simulation  Tools  for 
Virtual  Prototyping  of  Biomicrofluidic  Devices  and  Integrated  Systems”  focuses  on  the 
development  of  a  SPICE-like  environment  to  link  primitive  MEMS  elements  to  form  a 
system-level  model  for  coupled  fluid-thermal-chemical-electrical-magnetic-structural 
phenomena  in  microsystems.  This  project  will  validate  these  models  on  microfluidic  systems 
for  complex  chemical  and  biological  applications. 

CFDRC  will  continue  to  actively  commercialize  the  CFD-ACE+MEMS  software  to  the  MEMS 
community.  This  software  will  be  demonstrated  at  all  leading  MEMS  conferences,  workshops 
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and  trade-shows.  CFDRC  will  provide  academic  licenses  (at  a  highly  discounted  price)  to 
universities,  research  laboratories  and  non-profit  organizations  to  facilitate  widespread  use  of  the 
software.  Additionally,  CFDRC  will  also  provide  hands-on  training,  documentation,  technical 
support  and  upgrades  for  the  effective  use  of  the  software  (by  the  industry)  as  a  design  tool. 
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